Page 46 - Compact Numerical Methods For Computers

36 Compact numerical methods for computers
investigated by T Hoy Booker (Booker 1985). The proof in exact arithmetic is
incomplete. However, for a method such as the algorithm presented here, which uses
tolerances for zero, Booker has shown that the cyclic sweeps must eventually
terminate.

Algorithm 1. Singular-value decomposition

procedure 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 the
first 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 a
conventional 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. An
extended form of this commentary can be displayed on the screen by
removing the comment braces on the writeln statements below.

Copyright 1988 J. C. Nash
}
Var
i, 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 do
begin
D1 := W[i,j]; D2 := W[i,k];
W[i,j] := D1*c0+D2*s0; W[i,k] := -D1*s0+D2*c0
end; { 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’);