Nash - Compact Numerical Methods for Computers (523163), страница 13
Текст из файла (страница 13)
Solutions (a), (b), (c) and (d) in table 3.2 therefore present the solutionscorresponding to all, four, three and two principal components, respectively. Notethat these have 8, 9, 10 and 11 degrees of freedom because we estimate thecoefficients of the principal components, then transform these to give solutions interms of our original variables. The solution given by only three principalcomponents is almost as good as that for all components, that is.
a conventionalleast-squares solution. However, the coefficients in solutions (a), (b) and (c) arevery different.Neither the algorithms in this book nor those anywhere else can make a clearand final statement as to which solution is ‘best’. Here questions of statisticalsignificance will not be addressed, though they would probably enter into consideration if we were trying to identify and estimate a model intended for use inTABLE 3.2. Solutions for various principal-component regressions using the data in table 3.1.Tolerancefor zero(a)0(b)1(c)22(d)40R*0·972586(0·958879)0·969348(0·959131)0·959506(0·951407)0·93839(0·932789)X nitrogenxconstantx phosphatexcotashxpetroleum–0·290373–0·0461911·0194–0·159834·33368E–3–5·85314E–21·1757–0·2522960·6996215·14267E–34·34851E–20·392026–6·93389E–21·0115207·7822·54597E–3–0·152990·3001270·4692940·528881Regression involving a constant and index numbers for phosphate(e)0(f)10·968104(0·965204)0·273448179·3752·24851E–3The values in parentheses below each R’ are the corrected——0·5189660·945525statistic given by formula (3.52).—48Compact numerical methods for computerssome analysis or prediction.
To underline the difficulty of this task, merelyconsider the alternative modelincome = x1 + x3 (phosphate)(3.53)for which the singular values are computed as 1471·19 and 0·87188, again quitecollinear. The solutions are (e) and (f) in table 3.2 and the values of R2 speak forthemselves.A sample driver program DR0102.PAS is included on the program diskette.Appendix 4 describes the sample driver programs and supporting procedures andfunctions.Previous Home NextChapter 4HANDLING LARGER PROBLEMS4.1. INTRODUCTIONThe previous chapter used plane rotations multiplying a matrix from the right toorthogonalise its columns.
By the essential symmetry of the singular-value decomposition, there is nothing to stop us multiplying a matrix by plane rotations fromthe left to achieve an orthogonalisation of its rows. The amount of work involvedis of order m2 n operations per sweep compared to mn 2 for the columnwiseorthogonalisation (A is m by n), and as there are normally more rows thancolumns it may seem unprofitable to do this.
However, by a judicious combinationof row orthogonalisation with Givens’ reduction, an algorithm can be devisedwhich will handle a theoretically unlimited number of rows.4.2. THE GIVENS’ REDUCTIONThe above approach to the computation of a singular-value decomposition andleast-squares solution works very well on a small computer for relatively smallmatrices. For instance, it can handle least-squares regression calculations involving up to 15 economic time series with 25 years of data in less than 2000 words ofmain memory (where one matrix element occupies two words of storage). Whilesuch problems are fairly common in economics, biological and sociologicalphenomena are likely to have associated with them very large numbers ofobservations, m, even though the number of variables, n, may not be large.
Again,the formation of the sum-of-squares and cross-products matrix AT A and solutionvia the normal equations (2.22) should be avoided. In order to circumvent thedifficulties associated with the storage of the whole of the matrix A, the Givens’reduction can be used. A Givens’ transformation is simply a plane rotation whichtransforms two vectors so that an element of one of them becomes zero. Considertwo row vectors z T and yT and a pre-multiplying rotation:(4.1)wherec = cos φs = sin φ(4.2)and φ is the angle of rotation.
If Y1 is to be zero, then–sz1 + cy1 = 0(4.3)so that the angle of rotation in this case is given bytan φ = s/c = yl/z1.49(4.4)50Compact numerical methods for computersThis is a simpler angle calculation than that of §3.3 for the orthogonalisationprocess, since it involves only one square root per rotation instead of two. That is,if(4.5)then we haveandc = z 1 /p(4.6)s = y1/p.(4.7)It is possible, in fact, to perform such transformations with no square roots atall (Gentleman 1973, Hammarling 1974, Golub and Van Loan 1983) but no way hasso far come to light for incorporating similar ideas into the orthogonalising rotationof §3.3. Also, it now appears that the extra overhead required in avoiding the squareroot offsets the expected gain in efficiency, and early reports of gains in speed nowappear to be due principally to better coding practices in the square-root-freeprograms compared to their conventional counterparts.The Givens’ transformations are assembled in algorithm 3 to triangularise a realm by n matrix A.
Note that the ordering of the rotations is crucial, since anelement set to zero by one rotation must not be made non-zero by another.Several orderings are possible; algorithm 3 acts column by column, so thatrotations placing zeros in column k act on zeros in columns 1, 2, . .
. , (k - 1) andleave these elements unchanged. Algorithm 3 leaves the matrix A triangular, thatisA [i,j] = 0for i > j(4.8)which will be denoted R. The matrix Q contains the transformations, so that theoriginal m by n matrix isA = QR.(4.9)In words, this procedure simply zeros the last (m – 1) elements of column 1,then the last (m – 2) elements of column 2, .
. . , and finally the last ( m – n )elements of column n.Since the objective in considering the Givens’ reduction was to avoid storing alarge matrix, it may seem like a step backwards to discuss an algorithm whichintroduces an m by m matrix Q. However, this matrix is not needed for thesolution of least-squares problems except in its product Q T b with the right-handside vector b. Furthermore, the ordering of the rotations can be altered so thatthey act on one row at a time, requiring only storage for this one row and for theresulting triangular n by n matrix which will again be denoted R, that is(4.10)In the context of this decomposition, the normal equations (2.22) become(4.11)Handling larger problems51Thus the zeros below R multiply the last (m – n) elements of(4.12)where d1 is of order n and d2 of order (m – n).
ThusAT Ax = RT Rx= R T dl + 0d2 = AT b.(4.13)These equations are satisfied regardless of the values in the vector d2 bysolutions x to the triangular system(4.14)Rx = d1.This system is trivial to solve if there are no zero elements on the diagonal of R.Such zero elements imply that the columns of the original matrix are not linearlyindependent. Even if no zero or ‘small’ element appears on the diagonal, theoriginal data may be linearly dependent and the solution x to (4.14) in some way‘unstable’.Algorithm 3. Givens’ reduction of a real rectangular matrixprocedure givens(nRow,nCol : integer; (size of matrix to be decomposed)var A, Q: rmatrix); (the matrix to be decomposedwhich with Q holds the resulting decomposition){alg03.pas ==Givens’ reduction of a real rectangular matrix to Q * R form.This is a conventional version, which works one column at a time.Copyright 1988 J.
C. Nash}vari, j, k, mn: integer; {loop counters and array indices}b, c, eps, p, s : real; {angle parameters in rotations}beginwriteln(‘alg03.pas -- Givens’,chr(39),’ reduction -- column-wise’);{STEP 0 -- partly in the procedure call}mn := nRow; if nRow>nCol then mn := nCo1; {mn is the minimum of nRowand nCo1 and gives the maximum size of the triangular matrixresulting from the reduction. Note that the decomposition isstill valid when nRow<nCol, but the R matrix is thentrapezoidal, i.e.
an upper trianglular matrix with addedcolumns on the right.}for i := 1 to nRow dobegin {set Q to a unit matrix of size nRow}for j := 1 to nRow do Q[i,j] := 0.0;Q[i,i] := 1.0;end; {loop on i -- Q now a unit matrix}eps := calceps; {the machine precision}for j := 1 to (mn-1) do {main loop on diagonals of triangle} {STEP 1}begin {STEP 2}for k := (j+1) to nRow do {loop on column elements}begin {STEP 3}c := Au[j,j]; s := A[k,j]; {the two elements in the rotation}52Compact numerical methods for computersAlgorithm 3. Givens’ reduction of a real rectangular matrix (cont.)b := abs(c); if abs(s)>b then b := abs(s);if b>0 thenbegin {rotation is needed}c := c/b; s := s/b; {normalise elements to avoid over- or under-flow}p := sqrt(c*c+s*s); {STEP 4}s := s/p;if abs(s)>=eps then {STEP 5}begin {need to carry out rotations} {STEP 6}c := c/p;for i := 1 to nCo1 do {STEP 7 -- rotation of A}beginp := A[j,i]; A[j,i] := c*p+s*A[k,i]; A[k,i] := -s*p+c*A[k,i];end; {loop on i for rotation of A}for i := 1 to nRow do {STEP 8 -- rotation of Q.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
