Nash - Compact Numerical Methods for Computers (523163), страница 16
Текст из файла (страница 16)
pas -- Givens' reduction, svd, and least squares solutionOrder of 1s problem and no. of right hand sides = 51Enter a number to indicate end of data -9.9900000000E+02Obsn 1 563.00000 262.00000 461.00000 221.000001.00000 305.00000Uncorrelated residual(s) : 0.0000E+00Obsn 2 658.00000 291.00000 473.00000 222.000001.00000 342.00000Uncorrelated residual(s) : 0.0000E+00Obsn 3 676.00000 294.00000 513.00000 221.000001.00000 331.00000Uncorrelated residual(s) : 0.0000E+00Obsn 4 749.00000 302.00000 516.00000 218.000001.00000 339.00000Uncorrelated residual(s) : 0.0000E+00Obsn 5 834.00000 320.00000 540.00000 217.000001.00000 354.00000Uncorrelated residual(s) : 0.0000E+00Obsn 6 973.00000 350.00000 596.00000 218.000001.00000 369.00000Uncorrelated residual(s) : -6.563E-02Obsn 7 1079.00000 386.00000 650.00000 218.000001.00000 378.00000Uncorrelated residual(s) : -9.733E+00Obsn 8 1151.00000 401.00000 676.00000 225.000001.00000 368.00000Uncorrelated residual(s) : -6.206E+00Obsn 9 1324.00000 446.00000 769.00000 228.000001.00000 405.00000Uncorrelated residual(s) : 1.7473E+01Obsn 10 1499.00000 492.00000 870.00000 230.000001.00000 438.00000Uncorrelated residual(s) : 1.5054E+O1Obsn 11 1690.00000 510.00000 907.00000 237.000001.00000 438.00000Uncorrelated residual(s) : 7.4959E+00Obsn 12 1735.00000 534.00000 932.00000 235.000001.00000 451.00000Uncorrelated residual(s) : 1.0754E+00Obsn 13 1778.00000 559.00000 956.00000 236.000001.00000 485.00000Uncorrelated residual(s) : 1.5580E+O1Sweep 1 10 rotations performedSweep 2 10 rotations performedSweep 3 2 rotations performedSweep 4 0 rotations performedSingular values and principal componentsHandling larger problems63Singular value [1]=5.2985598853E+030.82043 0.27690 0.47815 0.14692 0.00065Singular value [2]=3.4551146213Et02-0.49538 0.30886 0.46707 0.66411 0.00322Singular value [3]=3.6112521703E+01-0.26021 -0.12171 0.71337 -0.63919 -0.00344Singular value [4]=2.1420869565E+010.11739 -0.90173 0.21052 0.35886 0.00093Singular value [5]=5.1382810120E-020.00006 -0.00075 0.00045 -0.00476 0.99999Enter a tolerance for zero (<0 to exit) 0.0000000000E+00Solution component [1]= -4.6392433678E-02Solution component [2]= 1.01938655593+00Solution component [3]= -1.5982291948E-01Solution component [4]= -2.9037627732E-01Solution component [5]= 2.0778262574Et02Residual sum of squares =9.6524564856E+02Enter a tolerance for zero (<0 to exit) 1.0000000000E+00Solution component [1]= -5.8532203918E-02Solution component [2]= 1.1756920631E+00Solution component [3]= -2.5228971048E-01Solution component [4]= 6.9962158969E-01Solution component [5]= 4.3336659982E-03Residual sum of squares =1.0792302647E+03Enter a tolerance for zero (<0 to exit) -1.0000000000E+004.5.
RELATED CALCULATIONSIt sometimes happens that a least-squares solution has to be updated as new dataare collected or become available. It is preferable to achieve this by means of astable method such as the singular-value decomposition. Chambers (1971) discusses the general problem of updating regression solutions, while Businger (1970)has proposed a method for updating a singular-value decomposition. However,the idea suggested in the opening paragraph of this chapter, in particular toorthogonalise (n + 1) rows each of n elements by means of plane rotations, worksquite well.
Moreover, it can be incorporated quite easily into algorithm 4, thougha little caution is needed to ensure the correct adjustment of quantities needed tocompute statistics such as R2. Nash and Lefkovitch (1977) present both FORTRANand BASIC programs which do this. These programs are sub-optimal in the sensethat they perform the normal sweep strategy through the rows of W, whereaswhen a new observation is appended the first n rows are already mutuallyorthogonal. Because the saving only applies during the first sweep, no specialsteps have been taken to employ this knowledge.
Unfortunately, each neworthogonalisation of the rows may take as long as the first, that is, the one thatfollows the Givens’ reduction. Perhaps this is not surprising since new observations may profoundly change the nature of a least-squares problem.64Compact numerical methods for computersThe method suggested is mainly useful for adding single observations, and otherapproaches are better if more than a very few observations are to be included. Forinstance, one could update the triangular form which results from the Givens’reduction if this had been saved, then proceed to the singular-value decompositionas in algorithm 4.No methods will be discussed for removing observations, since while methodsexist to accomplish this (see Lawson and Hanson 1974, pp 225-31), the operation is potentially unstable.
See also Bunch and Nielsen (1978).For instance, suppose we have a Givens’ QR decomposition of a matrix A (orany other QR decomposition with Q orthogonal and R upper-triangular), then addand delete a row (observation) denoted y T. Then after the addition of this row,the (1, 1) elements of the matrices are related by(4.31)where the tilde is used to indicate matrices which have been updated. Deletion ofy T now requires the subtraction(4.32)to be performed in some way or another, an operation which will involve digitcancellation if y1 andare close in magnitude.
The same difficulty may ofcourse occur in other columns-the first is simply easier to illustrate. Such casesimply that an element of yT dominates the column in which it occurs and as suchshould arouse suspicions about the data. Chambers’ (1971) subroutine to deleterows from a QR decomposition contains a check designed to catch such occurrences.Of interest to those users performing regression calculations are the estimates ofstandard errors of the regression coefficients (the least-squares solution elements).The traditional standard error formula isSE(bi ) = ( σ2(AT A) i i - 1 )½(4.33)where σ 2 is an estimate of the variance of data about the fitted model calculated bydividing the sum of squared residuals by the number of degrees of freedom (nRow –nCol) = (nRow – n).
The sum of squared residuals has already been computed inalgorithm 4, and has been adjusted for rank deficiency within the solution phase ofthe code.The diagonal elements of the inverse of the sum of squares and cross-productsmatrix may seem to pose a bigger task. However, the singular-value decompositionleads easily to the expressionT(ATA)-1 = VS + S+ V .(4.34)In particular, diagonal elements of the inverse of the sum of squares and cross-Handling larger problems65products matrix are(4.35)Thus, the relevant information for the standard errors is obtained by quite simplerow sums over the V matrix from a singular-value decomposition.
When the originalA matrix is rank deficient, and we decide (via the tolerance for zero used to select‘non-zero’ singular values) that the rank is r, the summation above reduces to(4.36)However, the meaning of a standard error in the rank-deficient case requires carefulconsideration, since the standard error will increase very sharply as small singularvalues are included in the summation given in (4.36). I usually refer to the dispersionmeasures computed via equations (4.33) through (4.36) for rank r < n cases as‘standard errors under the condition that the rank is5 (or whatever value r currentlyhas)‘.
More discussion of these issues is presented in Searle (1971) under the topic‘estimable functions’, and in various sections of Belsley, Kuh and Welsch (1980).Previous Home NextChapter 5SOME COMMENTS ON THE FORMATION OF THETCROSS-PRODUCTS MATRIX A ACommonly in statistical computations the diagonal elements of the matrix(ATA) -1(5.1)are required, since they are central to the calculation of variances for parametersestimated by least-squares regression. The cross-products matrix ATA from thesingular-value decomposition (2.53) is given byATA = VSUT USVT = VS2 V T .(5.2)TThis is a singular-value decomposition of A A , so that(ATA)+ = V(S+ )2 V T .(5.3)If the cross-products matrix is of full rank, the generalised inverse is identical tothe inverse (5.1) and, further,S+ = S- 1 .(5.4)Thus we have(AT A)-1 = VS- 2VT .(5.5)The diagonal elements of this inverse are therefore computed as simple rownorms of the matrixVS- l .(5.6)In the above manner the singular-value decomposition can be used to computethe required elements of the inverse of the cross-products matrix.
This means thatthe explicit computation of the cross-products matrix is unnecessary.Indeed there are two basic problems with computation of AT A. One is inducedby sloppy programming practice, the other is inherent in the formation of AT A.The former of these occurs in any problem where one of the columns of A isconstant and the mean of each column is not subtracted from its elements. Forinstance, let one of the columns of A (let it be the last) have all its elements equalto 1.
The normal equations (2.22) then yield a cross-products matrix with last row(and column), say the nth,(5.7)But(5.8)66Some comments on the formation of the cross-products matrix67where is the mean of the jth column of the m by n matrix A. Furthermore, theright-hand side of the nth normal equation is(5.9)This permits xn to be eliminated by using the nth normal equation(5.10)or(5.11)When this expression is substituted into the normal equations, the kth equation(note carefully the bars above the symbols) becomes(5.12)But since(5.13)and(5.14)equation (5.12) becomes(5.15)which defines a set of normal equations of order (n - 1)(A' ) TA' x' = (A') Tb'(5.16)where A' is formed from the (n – 1) non-constant columns of A each adjusted bysubtraction of the mean and where b' is formed from b by subtraction of .
x' issimply x without xn.Besides reducing the order of the problem, less information is lost in theformation of (A') TA' than AT A, since the possible addition of large numbers to thematrix is avoided. These large numbers have subsequently to be subtracted fromeach other in the solution process, and this subtraction leads to digit cancellationwhich one should always seek to avoid.As an example, consider the calculation of the variance of the column vector(5.17)68Compact numerical methods for computerswhich has mean= 1003·5 so that(5.18)The variance is computed via either(5.19)where m = 4 is the number of elements in the vector, or since(5.20)by(5.21)Note that statisticians prefer to divide by (m – 1) which makes a change necessaryin (5.21). Equation (5.19) when applied to the example on a six decimal digitcomputer givesvar(a) = (12·25 + 2·25 + 0·25 + 20·25)/4 = 35/4 = 8·75.(5.22)By comparison, formula (5.21) produces the results in table 5.1 depending onwhether the computer truncates (chops) or rounds.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
