Page 47 - Compact Numerical Methods For Computers

Singular-value decomposition, and use in least-squares problems 37
Algorithm 1. Singular-value decomposition (cont.)
writeln(‘of Z. The final array W contains the decomposition in a’);
writeln(‘special form, namely,’);
writeln;
writeln(‘ ( U S ) ’);
writeln(‘ W = ( ) ’);
writeln(‘ ( V ) ’);
writeln;
writeln(‘The matrices U and V are extracted from W, and S is’);
writeln(‘found from Z. However, the (U S) matrix and V matrix may’);
writeln(‘also be used directly in calculations, which we prefer’);
writeln(‘since fewer arithmetic operations are then needed.’);
writeln;

{STEP 0 Enter nRow, nCo1, the dimensions of the matrix to be decomposed.}
writeln(‘alg01.pas--NashSVD’);
eps := Calceps; {Set eps, the machine precision.}
slimit := nCo1 div 4; if slimit< then slimit := 6;
{Set slimit, a limit on the number of sweeps allowed. A suggested
limit is max([nCol/4], 6).}
SweepCount := 0; {to count the number of sweeps carried out}
e2 := 10.0*nRow*eps*eps;
tol := eps*0.1;
{Set the tolerances used to decide if the algorithm has converged.
For further discussion of this, see the commentary under STEP 7.}
EstColRank := nCo1; {current estimate of rank};
{Set V matrix to the unit matrix of order nCo1.
V is stored in rows (nRow+1) to (nRow+nCol) of array W.}
for i := 1 to nCo1 do
begin
for j := 1 to nCo1 do
W[nRow+i,j] := 0.0; W[nRow+i,i] := 1.0;
end; {loop on i, and initialization of V matrix}
{Main SVD calculations}
repeat {until convergence is achieved or too many sweeps are carried out}
RotCount := EstColRank*(EstColRank-1) div 2; {STEP 1 -- rotation counter}
SweepCount := SweepCount+1;
for j := 1 to EstColRank-1 do {STEP 2 -- main cyclic Jacobi sweep}
begin {STEP 3}
for k := j+l to EstColRank do
begin {STEP 4}
p := 0.0; q := 0.0; r := 0.0;
for i := 1 to nRow do {STEP 5}
begin
x0 := W[i,j]; y0 := W[i,k];
p := p+x0*y0; q := q+x0*x0; r := r+y0*y0;
end;
Z[j] := q; Z[k] := r;
{Now come important convergence test considerations. First we
will decide if rotation will exchange order of columns.}
if q >= r then {STEP 6 -- check if the columns are ordered.}
begin {STEP 7 Columns are ordered, so try convergence test.}
if (q<=e2*2[1]) or (abs(p)<= tol*q) then RotCount := RotCount-1
{There is no more work on this particular pair of columns in the