Nash - Compact Numerical Methods for Computers (523163), страница 11
Текст из файла (страница 11)
If x = ak, then straightforward use of equations (3.23) and (3.24)or (3.25) and (3.26) givesXT X – x T x = (v – q)/2 > 0.(3.33)Using (3.31) and (3.22) in (3.33) gives(3.34)Thus, once columns become sufficiently separated by size and the nonorthogonality sufficiently diminished, the column ordering is stable. When somecolumns are equal in norm but orthogonal, the above theorem can be applied tocolumns separated by size.The general question of convergence in the case of equal singular values has been36Compact numerical methods for computersinvestigated by T Hoy Booker (Booker 1985). The proof in exact arithmetic isincomplete.
However, for a method such as the algorithm presented here, which usestolerances for zero, Booker has shown that the cyclic sweeps must eventuallyterminate.Algorithm 1. Singular-value decompositionprocedure NashSVD(nRow, nCo1: integer; {size of problem}var W: wmatrix; {working matrix}var Z: rvector); {squares of singular values}(alg01.pas ==form a singular value decomposition of matrix A which is stored in thefirst nRow rows of working array W and the nCo1 columns of this array.The first nRow rows of W will become the product U * S of aconventional svd, where S is the diagonal matrix of singular values.The last nCo1 rows of W will be the matrix V of a conventional svd.On return, Z will contain the squares of the singular values.
Anextended form of this commentary can be displayed on the screen byremoving the comment braces on the writeln statements below.Copyright 1988 J. C. Nash}Vari, j, k, EstColRank, RotCount, SweepCount, slimit : integer;eps, e2, tol, vt, p, h2, x0, y0, q, r, c0, s0, c2, d1, d2 : real;procedure rotate; (STEP 10 as a procedure}(This rotation acts on both U and V, by storing V at the bottom of U}begin (<< rotation )for i := 1 to nRow+nCol dobeginD1 := W[i,j]; D2 := W[i,k];W[i,j] := D1*c0+D2*s0; W[i,k] := -D1*s0+D2*c0end; { rotation >>}end; { rotate }begin { procedure SVD }{ -- remove the comment braces to allow message to be displayed -writeln(‘Nash Singular Value Decomposition (NashSVD).’);writeln;writeln(‘The program takes as input a real matrix A.’);writeln;writeln(‘Let U and V be orthogonal matrices, & S’);writeln(‘a diagonal matrix, such that U” A V = S.’);writeln(‘Then A = U S V” is the decomposition.’);writeln(‘A is assumed to have more rows than columns.
If it’);writeln(‘does not, the svd of the transpose A” gives the svd’);writeln(‘of A, since A” = V S U”.’);writeln;writeln(‘If A has nRow rows and nCo1 columns, then the matrix’);writeln(‘is supplied to the program in a working array W large’);writeln(‘enough to hold nRow+nCol rows and nCo1 columns.’);writeln(‘Output comprises the elements of Z, which are the’);writeln(‘squares of the elements of the vector S, together’);writeln(‘with columns of W that correspond to non-zero elements’);Singular-value decomposition, and use in least-squares problemsAlgorithm 1.
Singular-value decomposition (cont.)writeln(‘of Z. The final array W contains the decomposition in a’);writeln(‘special form, namely,’);writeln;writeln(‘ ( U S ) ’);writeln(‘ W = ( ) ’);writeln(‘ ( V ) ’);writeln;writeln(‘The matrices U and V are extracted from W, and S is’);writeln(‘found from Z. However, the (U S) matrix and V matrix may’);writeln(‘also be used directly in calculations, which we prefer’);writeln(‘since fewer arithmetic operations are then needed.’);writeln;{STEP 0 Enter nRow, nCo1, the dimensions of the matrix to be decomposed.}writeln(‘alg01.pas--NashSVD’);eps := Calceps; {Set eps, the machine precision.}slimit := nCo1 div 4; if slimit< then slimit := 6;{Set slimit, a limit on the number of sweeps allowed.
A suggestedlimit is max([nCol/4], 6).}SweepCount := 0; {to count the number of sweeps carried out}e2 := 10.0*nRow*eps*eps;tol := eps*0.1;{Set the tolerances used to decide if the algorithm has converged.For further discussion of this, see the commentary under STEP 7.}EstColRank := nCo1; {current estimate of rank};{Set V matrix to the unit matrix of order nCo1.V is stored in rows (nRow+1) to (nRow+nCol) of array W.}for i := 1 to nCo1 dobeginfor j := 1 to nCo1 doW[nRow+i,j] := 0.0; W[nRow+i,i] := 1.0;end; {loop on i, and initialization of V matrix}{Main SVD calculations}repeat {until convergence is achieved or too many sweeps are carried out}RotCount := EstColRank*(EstColRank-1) div 2; {STEP 1 -- rotation counter}SweepCount := SweepCount+1;for j := 1 to EstColRank-1 do {STEP 2 -- main cyclic Jacobi sweep}begin {STEP 3}for k := j+l to EstColRank dobegin {STEP 4}p := 0.0; q := 0.0; r := 0.0;for i := 1 to nRow do {STEP 5}beginx0 := W[i,j]; y0 := W[i,k];p := p+x0*y0; q := q+x0*x0; r := r+y0*y0;end;Z[j] := q; Z[k] := r;{Now come important convergence test considerations.
First wewill decide if rotation will exchange order of columns.}if q >= r then {STEP 6 -- check if the columns are ordered.}begin {STEP 7 Columns are ordered, so try convergence test.}if (q<=e2*2[1]) or (abs(p)<= tol*q) then RotCount := RotCount-1{There is no more work on this particular pair of columns in the3738Compact numerical methods for computersAlgorithm 1. Singular-value decomposition (cont.)current sweep. That is, we now go to STEP 11. The firstcondition checks for very small column norms in BOTH columns, forwhich no rotation makes sense.
The second condition determinesif the inner product is small with respect to the larger of thecolumns, which implies a very small rotation angle.)else {columns are in order, but their inner product is not small}begin {STEP 8}p := p/q; r := 1-r/q; vt := sqrt(4*p*p + r*r);c0 := sqrt(0.5*(1+r/vt)); s0 := p/(vt*c0);rotate;endend {columns in order with q>=r}else {columns out of order -- must rotate}begin {STEP 9}{note: r > q, and cannot be zero since both are sums of squares forthe svd. In the case of a real symmetric matrix, this assumptionmust be questioned.}p := p/r; q := q/r-1; vt := sqrt(4*p*p + q*q);s0 := sqrt(0.5*(1-q/vt));if p<0 then s0 := -s0;co := p/(vt*s0);rotate; {The rotation is STEP 10.}end;{Both angle calculations have been set up so that large numbers donot occur in intermediate quantities.
This is easy in the svd case,since quantities x2,y2 cannot be negative. An obvious scaling forthe eigenvalue problem does not immediately suggest itself.)end; {loop on K -- end-loop is STEP 11}end; {loop on j -- end-loop is STEP 12}writeln(‘End of Sweep #’, SweepCount,‘- no.
of rotations performed =’, RotCount);{STEP 13 -- Set EstColRank to largest column index for whichZ[column index] > (Z[1]*tol + tol*tol)Note how Pascal expresses this more precisely.}while (EstColRank >= 3) and (Z[EstColRank] <= Z[1]*tol + tol*tol)do EstColRank := EstColRank-1;{STEP 14 -- Goto STEP 1 to repeat sweep if rotations have beenperformed and the sweep limit has not been reached.}until (RotCount=0) or (SweepCount>slimit);{STEP 15 -- end SVD calculations}if (SweepCount > slimit) then writeln(‘**** SWEEP LIMIT EXCEEDED’);if (SweepCount > slimit) then{Note: the decomposition may still be useful, even if the sweeplimit has been reached.}end; {alg01.pas == NashSVD}3.5.
AN ALTERNATIVE IMPLEMENTATION OF THE SINGULARVALUE DECOMPOSITIONOne of the most time-consuming steps in algorithm 1 is the loop which comprisesSingular-value decomposition, and use in least-squares problems39TSTEP 5. While the inner product used to compute p = x y must still be performed,it is possible to use equation (3.33) and the corresponding result for Y, that isY T Y – yT y = – ( v – q)/2(3.35)to compute the updated column norms after each rotation.
There is a danger thatnearly equal magnitudes may be subtracted, with the resultant column norm having alarge relative error. However, if the application requires information from the largestsingular values and vectors, this approach offers some saving of effort. The changesneeded are:(1) an initial loop to compute the Z[i], that is, the sum of squares of the elementsof each column of the original matrix A;(2) the addition of two statements to the end of the main svd loop on k, which,if a rotation has been performed, update the column norms Z[j] and Z[k] viaformulae (3.34) and (3.35).
Note that in the present algorithm the quantitiesneeded for these calculations have not been preserved. Alternatively, add atthe end of STEP 8 (after the rotation) the statementsZ[j] : = Z[j] + 0·5*q*(vt – r);Z[k] : = Z[k] – 0·5*q*(vt – r);if Z[k] < 0·0 then Z[k] : = 0·0;and at the end of STEP 9 the statementsZ[j] := Z[j] + 0·5*r*(vt – q);Z[k] := Z[k] – 0·5*r*(vt – q);if Z[k] < 0·0 then Z[k] := 0·0;(3) the deletion of the assignmentsZ[j] := q;Z[k] := r;at the end of STEP 5.As an illustration of the consequences of such changes, the singular-value decompositions of a 6 by 4 matrix derived from the Frank symmetric test matrix and an8 by 5 portion of the Hilbert matrix were calculated. In the latter test the Turbo-87Pascal compiler was used rather than the regular Turbo Pascal compiler (versions3.01a of both systems).
The results below present the modified algorithm result(s)above the corresponding figure(s) for the regular method.Frank Matrix:Column orthogonality of ULargest inner product is 4, 4 = -4.7760764232E-09Largest inner product is 4, 4 = 3.4106051316E-1240Compact numerical methods for computersSingular values3.3658407311E+00 1.0812763036E+00 6.7431328720E-01 5.3627598567E-013.3658407311E+00 1.0812763036E+00 6.7431328701E-01 5.3627598503E-01Hilbert segment:Column orthogonality of ULargest inner product is 5, 5 = -1.44016460160157E-006Largest inner product is 3, 3 = 5.27355936696949E-016Singular values1.27515004411E+000 4.97081651063E-001 1.30419686491E-001 2.55816892287E-0021.27515004411E+000 4.97081651063E-001 1.30419686491E-001 2.55816892259E-0023.60194233367E-0033.60194103682E-0033.6.
USING THE SINGULAR-VALUE DECOMPOSITION TO SOLVELEAST-SQUARES PROBLEMSBy combining equations (2.33) and (2.56), the singular-value decomposition canbe used to solve least-squares problems (2.14) viax = VS+ UT b.(3.36)+However, the definition (2.57) of S is too strict for practical computation,since a real-world calculation will seldom give singular values which are identically zero. Therefore, for the purposes of an algorithm it is appropriate to define(3.37)where q is some tolerance set by the user.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
