Nash - Compact Numerical Methods for Computers (523163), страница 14
Текст из файла (страница 14)
Note: nRow not nCo1.}beginp := Q[i,j]; Q[i,jl := c*p+s*Q[i,k]; Q[i,k] := -s*p+c*Q[i,k];end; {loop on i for rotation of Q}end; {if abs(s)>=eps}end; {if b>0}end; {loop on k -- the end-loop is STEP 9}end; {loop on j -- the end-loop is STEP 10}end; {alg03.pas == Givens’ reduction}After the Givens’ procedure is complete, the array A contains the triangular factor R in rows 1to mn. If nRow is less than nCol, then the right-hand (nCo1 – nRow) columns of array A contain thetransformed columns of the original matrix so that the product Q*R = A, in which R is nowtrapezoidal. The decomposition can be used together with a back-substitution algorithm such asalgorithm 6 to solve systems of linear equations.The order in which non-zero elements of the working array are transformed to zero is notunique.
In particular, it may be important in some applications to zero elements row by rowinstead of column by column. The ,file alg03a.pas on the software disk presents such a row-wisevariant. Appendix 4 documents driver programs DR03.PAS and DR03A.PAS which illustratehow the two Givens’ reduction procedures may be used.Example 4.1. The operation of Givens’ reductionThe following output of a Data General ECLIPSE operating in six hexadecimaldigit arithmetic shows the effect of Givens’ reduction on a rectangular matrix. Ateach stage of the loops of steps 1 and 2 of algorithm 3 the entire Q and A matricesare printed so that the changes are easily seen.
The loop parameters j and k aswell as the matrix elements c = A[j,j] and s = A[k,j] are printed also. In thisexample, the normalisation at step 3 of the reduction is not necessary, and thesine and cosine of the angle of rotation could have been determined directly fromthe formulae (4.5), (4.6) and (4.7).The matrix chosen for this example has only rank 2. Thus the last row of theFINAL A MATRIX is essentially null. In fact, small diagonal elements of thetriangular matrix R will imply that the matrix A is ‘nearly’ rank-deficient.However, the absence of small diagonal elements in R, that is, in the final array A,do not indicate that the original A is of full rank. Note that the recombination of53Handling larger problemsthe factors Q and R gives back the original matrix apart from very small errorswhich are of the order of the machine precision multiplied by the magnitude ofthe elements in question.*RUNTEST GIVENS - GIFT - ALG 3 DEC: 12 77SIZE -- M= ? 3 N= ? 4MTIN - INPUT M BY N MATRIXROW 1 : ? 1 ? 2 ? 3 ? 4ROW 2 : ? 5 ? 6 ? 7 ? 8ROW 3 : ? 9 ? 10 ? 11 ? 12ORIGINAL A MATRIX1ROW 1 :ROW 2 :59HOW 3 :26103711GIVENS TRIANGULARIZATION DEC: 12 77Q MATRIX10ROW 1 :01ROW 2 :0ROW 3 :0J= 1K= 2A MATRIXROW 1 :ROW 2 :ROW 3 :QROWROWROWMATRIX1 :2 :3 :K= 3J = 1A MATRIXROW 1 :ROW 2 :ROW 3 :QROWROWROWMATRIX1 :2 :3 :J= 2FINALROW 1ROW 2ROW 3A[J,J]= 15.09902-1.19209E-079.196116.9805810001A[K,J]= 56.27572-.78446610-.980581.l961160A[J,J]= 5.099027.45242-1.56893118.62912-2.353412001A[K,J]= 910.3441-1.19209E-07011.7942-.784466-.53086213.2443-1.56893-1.061729.66738E-02.483369.870063-.980581.1961160-.170634-.853168.492941K= 3A[J,J]=-.784466A[K,J]=-.530862A MATRIX11.794213.244310.3441:.947208:9.87278E-081.89441-6.68l09E-08:0-9.53674E-07FINAL Q MATRIXROW 1 :9.66738E-02ROW 2 :.483369ROW 3 :.870063.907738.315737.276269-.40825.816498-.4408249RECOMBINATIONROW 1 :ROW 2 :ROW 3 :2.000016.00002103.000017.000021115.000019.00001481214.6944-2.3534-1.5925814.69442.84162-1.90735E-064.000018.000021254Compact numerical methods for computers4.3.
EXTENSION TO A SINGULAR-VALUE DECOMPOSITIONIn order to determine if linear dependencies are present, it is possible to extendthe Givens’ reduction and compute a singular-value decomposition by orthogonalising the rows of R by plane rotations which act from the left. It is notnecessary to accumulate the transformations in a matrix U; instead they can beapplied directly to Q T b, in fact, to the first n elements of this vector which formd1.
The rotations can be thought of as acting all at once as the orthogonal matrixP T. Applying this to equation (4.14) givesPT R x = PT d1 = f.(4.15)THowever, the rows of P R are orthogonal, that isPT R = SVT(4.16)withSVT VS = S2 .(3.17)Combining the Givens’ reduction and orthogonalisation steps givesP T Q T A = PT(4.18)or(4.19)which is a singular-value decomposition.As in the case of the columnwise orthogonalisation, small singular values (i.e.rows of P TR having small norm) will cause V to possess some unnormalised rowshaving essentially zero elements. In this case (4.17) will not be correct. since(4.20)where k is the number of singular values larger than some pre-assigned tolerancefor zero.
Since in the solution of least-squares problems these rows always actonly in products with S or S+, this presents no great difficulty to programming analgorithm using the above Givens’ reduction/row orthogonalisation method.4.4. SOME LABOUR-SAVING DEVICESThe above method is not nearly so complicated to implement as it may appear.Firstly, all the plane rotations are row-wise for both the Givens’ reduction and theorthogonalisation. Moreover, one or more (say g) vectors b can be concatenatedwith the matrix A so that rotations do not have to be applied separately to these,but appear to act on a single matrix.The second observation which reduces the programming effort is that the rowsof this matrix (A, b) are needed only one at a time.
Consider a working arrayHandling larger problems55(n + 1) by (n + g) which initially has all elements in the first n rows equal to zero.Each of the m observations or rows of (A, b) can be loaded in succession into the(n + 1)th row of the working array. The Givens’ rotations will suitably fill up theworkspace and create the triangular matrix R. The same workspace suffices for theorthogonalisation.
Note that the elements of the vectors d2 are automatically leftin the last g elements of row (n + 1) of the working array when the first n havebeen reduced to zero. Since there are only (m – n) components in each d2, but mrows to process, at least n of the values left in these positions will be zero. Thesewill not necessarily be the first n values.A further feature of this method is that the residual sum of squares, r T r, isThis can be shown quite easily sinceequal to the sum of squared termsr T r = (b – Ax) T (b – Ax)= bTb – bT Ax – xTA Tb + xTATAx.(4.21)By using the normal equations (2.22) the last two terms of this expression cancelleaving(4.22)r T r = b T b – b T Ax.If least-squares problems with large numbers of observations are being solvedvia the normal equations, expression (4.22) is commonly used to compute theresidual sum of squares by accumulating bT b, AT A and AT b with a single passthrough the data.
In this case, however, (4.22) almost always involves thesubtraction of nearly equal numbers. For instance, when it is possible to approximate b very closely with Ax, then nearly all the digits in b T b will be cancelled bythose in b TAx, leaving a value for r T r with very few correct digits.For the method using rotations, on the other hand, we have(4.23)and(4.24)by equation (4.12).
Hence, by substitution of (4.23) and (4.24) into (4.22) weobtain(4.25)The cancellation is now accomplished theoretically with the residual sum ofsquares computed as a sum of positive terms, avoiding the digit cancellation.The result (4.25) is derived on the assumption that (4.14) holds. In therank-deficient case, as shown by k zero or ‘small’ singular values, the vector f inequation (4.15) can be decomposed so that(4.26)where f1 is of order (n – k) and f 2 of order k. Now equation (4.24) will have theform(4.27)56Compact numerical methods for computersby application of equation (4.16) and the condition that Sk+l, Sk+2, . .
. , Sn, are all‘zero’. Thus, using(4.28)and (4.22) with (4.27) and (4.23), the residual sum of squares in the rank-deficientcase is(4.29)From a practical point of view (4.29) is very convenient, since the computation ofthe residual sum of squares is now clearly linked to those singular values whichare chosen to be effectively zero by the user of the method. The calculation isonce again as a sum of squared terms, so there are no difficulties of digitcancellation.The vector(4.30)in the context of statistical calculations is referred to as a set of uncorrelatedresiduals (Golub and Styan 1973).Nash and Lefkovitch (1976) report other experience with algorithm 4. Inparticular, the largest problem I have solved using it involved 25 independentvariables (including a constant) and two dependent variables, for which there were196 observations.
This problem was initially run on a Hewlett-Packard 9830calculator, where approximately four hours elapsed in the computation. Later thesame data were presented to a FORTRAN version of the algorithm on both Univac1108 and IBM 370/168 equipment, which each required about 12 seconds ofprocessor time. Despite the order-of-magnitude differences in the timings between the computers and the calculator, they are in fact roughly proportional tothe cycle times of the machines. Moreover, as the problem has a singularity oforder 2, conventional least-squares regression programs were unable to solve it,and when it was first brought to me the program on the HP 9830 was the only oneon hand which could handle it.Algorithm 4.
Givens’ reductions, singular-value decomposition and least-squaressolutionprocedure GivSVD( n : integer; {order of problem}nRHS: integer; {number of right hand sides}var B: x-matrix; {matrix of solution vectors}var rss: r-vector; {residual sums of squares}var svs: x-vector; {singular values}var W: rmatrix; {returns V-transpose}var nobs : integer); {number of observations}{alg04.pas ==Givens’ reduction, singular value decomposition and least squaressolution.In this program, which is designed to use a very small working array yetsolve least squares problems with large numbers of observations, we do notexplicitly calculate the U matrix of the singular value decomposition.Handling larger problems57Algorithm 4. Givens’ reductions, singular-value decomposition and least-squares solution (cont.)One could save the rotations and carefully combine them to produce the Umatrix.
However, this algorithm uses plane rotations not only to zeroelements in the data in the Givens’ reduction and to orthogonalize rows ofthe work array in the svd portion of the code, but also to move the datainto place from the (n+1)st row of the working array into which the datais read. These movements i.e. of the observation number nobs,would normally move the data to row number nobs of the originalmatrix A to be decomposed. However, it is possible, as in the array givenby data file ex04.cnm3 1 <--- there are 3 columns in the matrix Aand 1 right hand side-999 <--- end of data flag1 2 3 1 <--- the last column is the RHS vector247122215311-999 0 0 0 <--- end of data rowthat this movement does not take place.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
