Page 48 - Compact Numerical Methods For Computers
P. 48

38                Compact numerical methods for computers
                            Algorithm 1. Singular-value decomposition (cont.)

                                           current sweep. That is, we now go to STEP 11. The first
                                           condition checks for very small column norms in BOTH columns, for
                                           which no rotation makes sense. The second condition determines
                                           if the inner product is small with respect to the larger of the
                                           columns, 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;
                                           end
                                        end {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 for
                                           the svd. In the case of a real symmetric matrix, this assumption
                                           must 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 do
                                        not occur in intermediate quantities. This is easy in the svd case,
                                         since quantities x2,y2 cannot be negative. An obvious scaling for
                                         the 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 which
                                      Z[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 been
                                         performed 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 sweep
                                limit has been reached.}
                             end; {alg01.pas == NashSVD}




                             3.5. AN ALTERNATIVE IMPLEMENTATION OF THE SINGULAR-
                                                  VALUE DECOMPOSITION
                           One of the most time-consuming steps in algorithm 1 is the loop which comprises
   43   44   45   46   47   48   49   50   51   52   53