Page 138 - Compact Numerical Methods For Computers
P. 138

Real symmetric matrices                    127

                     that
                                                                                       (10.41)
                                               (k+1)                          ( k )
                     so that each rotation causes A  to be ‘more diagonal’ than A  .
                       Specification of the order in which the off-diagonal elements are to be made
                     zero defines a particular variant of the Jacobi algorithm. One of the simplest is to
                     take the elements in row-wise fashion: (1, 2), (1, 3), . . . , (1, n), (2, 3), (2, 4), . . . ,
                     (2, n), . . . , (n – 1, n ). Such cyclic Jacobi algorithms have not been proved to
                     converge in general, except in the case where the angle of rotation f is
                     constrained so that
                                                  - p/4< f < p/4.                      (10.42)

                       Similarly to the orthogonalisation algorithm of §3.3, the difficulty lies in
                     demonstrating that the ordering of the diagonal elements is stable. For f
                     restricted as in (10.42), Forsythe and Henrici (1960) have carried out the
                     complicated proof, and most authors (Wilkinson 1965, Rutishauser 1966,
                     Schwarz et al 1973, Ralston 1965) discuss algorithms using angle calculations
                     which satisfy the restriction (10.42). In fact, among the variety of texts available
                     on numerical linear algebra, the author has noted only one (Fröberg 1965) which
                     uses the calculation based on equations (10.38), (10.39) and (3.22)-(3.27). The
                     advantage of the calculation given here is that the eigenvalues, which are
                     produced as the diagonal elements of the matrix that is the limit of the sequence
                       (k)
                     A , are ordered from most positive to most negative. Most applications which
                     require eigensolutions also require them ordered in some way and the ordering
                     that the current process yields is the one I have almost invariably been required to
                     produce. No extra program code is therefore needed to order the eigenvalues and
                     eigenvectors and there may be some savings in execution time as well.
                        We have already noted in chapter 3 the research of Booker (1985) relating to the
                      convergence of the cyclic Jacobi method with the ordering angle calculation.
                        Jacobi (1846) did not use a cyclic pattern, but chose to zero the largest
                      off-diagonal element in the current matrix. This process has generally been
                      considered inappropriate for automatic computation due to the time required for
                      the search before each plane rotation. However, for comparison purposes I
                      modified a BASIC precursor of algorithm 14. The changes made were only as many as
                      were needed to obtain the zeroing of the largest off-diagonal element in the
                      present matrix, and no attempt was made to remove computations present in
                      the algorithm which were used to provide convergence test information. The
                      processor time required for the order-5 Moler matrix (appendix 1) on a Data
                      General NOVA operating in six hexadecimal digit arithmetic was 15·7 seconds
                      for algorithm 14 while the Jacobi method required 13·7 seconds. Indeed the latter
                      required only 30 rotations while algorithm 14 took 5 sweeps (up to 50 rotations).
                      The comparison of timings on this one example may be misleading, especially as
                      the system uses an interpreter, that is, programs are executed by translating the
                      BASIC  at run time instead of compiling it first. (This has an analogy in the
                      translation of a speech either in total or simultaneously.) However, for matrices
                      with only a few large off-diagonal elements, the relatively simple changes to
   133   134   135   136   137   138   139   140   141   142   143