Page 44 - Compact Numerical Methods For Computers

34 Compact numerical methods for computers

Because only two columns are involved in the kth rotation, we have
2
T
2
T
Z (k) = Z (k-1) + (X Y) – (x y) .
(3.18)
But condition (3.14) implies
(k) (k-1) T 2
Z = Z – (x y)
(3.19)
so that the non-orthogonality is reduced at each rotation.
The specific formulae for the sine and cosine of the angle of rotation are (see
e.g. Nash 1975) given in terms of the quantities
T
p = x y (3.20)
T
T
q = x x – y y (3.21)
and
2 2 ½
v = (4p + q ) . (3.22)
They are
½
cos f = [(v + q)/(2v ) ] (3.23)
for q > 0
sin f = p/(v cos f ) (3.24)
½
sin f = sgn(p)[(v – q)/(2v) ] (3.25)
for q < 0
cos f = p/(u sin f) (3.26)
where
sgn (p) =} 1 for p > 0 (3.27)
–1 for p < 0.
Note that having two forms for the calculation of the functions of the angle of
rotation permits the subtraction of nearly equal numbers to be avoided. As the
matrix nears orthogonality p will become small, so that q and v are bound to have
nearly equal magnitudes.
In the first edition of this book, I chose to perform the computed rotation only
when q > r, and to use

sin ( f ) = 1 cos ( f ) = 0 (3.28)

when q < 0. This effects an interchange of the columns of the current matrix A.
However, I now believe that it is more efficient to perform the rotations as defined in
the code presented. The rotations (3.28) were used to force nearly null columns of the
final working matrix to the right-hand side of the storage array. This will occur when
the original matrix A suffers from linear dependencies between the columns (that is,
is rank deficient). In such cases, the rightmost columns of the working matrix
eventually reflect the lack of information in the data in directions corresponding to
the null space of the matrix A. The current methods cannot do much about this lack
of information, and it is not sensible to continue computations on these columns. In
the current implementation of the method (Nash and Shlien 1987), we prefer to
ignore columns at the right of the working matrix which become smaller than a