Page 49 - Compact Numerical Methods For Computers
P. 49

Singular-value decomposition, and use in least-squares problems  39
                                                                      T
                      STEP 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 is


                                                 T     T
                                                Y Y – y y = – ( v  – q)/2                (3.35)

                      to compute the updated column norms after each rotation. There is a danger that
                      nearly equal magnitudes may be subtracted, with the resultant column norm having a
                      large relative error. However, if the application requires information from the largest
                      singular values and vectors, this approach offers some saving of effort. The changes
                       needed are:
                          (1) an initial loop to compute the Z[i], that is, the sum of squares of the elements
                              of 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] via
                              formulae (3.34) and (3.35). Note that in the present algorithm the quantities
                              needed for these calculations have not been preserved. Alternatively, add at
                              the end of STEP 8 (after the rotation) the statements

                                  Z[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 statements

                                  Z[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 assignments

                                  Z[j] := q;    Z[k] := r;

                              at the end of STEP 5.
                         As an illustration of the consequences of such changes, the singular-value de-
                       compositions of a 6 by 4 matrix derived from the Frank symmetric test matrix and an
                       8 by 5 portion of the Hilbert matrix were calculated. In the latter test the Turbo-87
                       Pascal compiler was used rather than the regular Turbo Pascal compiler (versions
                       3.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 U
                       Largest inner product is 4, 4 = -4.7760764232E-09
                       Largest inner product is 4, 4 =  3.4106051316E-12
   44   45   46   47   48   49   50   51   52   53   54