Page 53 - Compact Numerical Methods For Computers

Singular-value decomposition, and use in least-squares problems 43

Algorithm 2. Least-squares solution via singular-value decomposition (cont.)
begin
writeln(‘alg02.pas == svdlss’);
repeat
writeln;
writeln(‘Singular values’);
for j := 1 to nCo1 do
begin
write(sqrt(Z[j]):18,‘’);
if j = 4 * (j div 4) then writeln;
end;
writeln;
write(‘Enter a tolerance for zero singular value (<0 to quit)’);
readln(infile,q);
if length(infname)>0 then writeln(q);
if q>=0.0 then
begin
q := q*q; {we will work with the Z vector directly}
for i := 1 to nCo1 do {loop to generate each element of Bvec}
begin
s := 0.0;
for j := 1 to nCo1 do {loop over columns of V}
begin
for k := 1 to nRow do {loop over elements of Y}
begin
if Z[j]>q then
s := s + W[i+nRow,j]*W[k,j]*Y[k]/Z[j];
{this is V * S+ * U-transpose * Y = A+ * Y}
{NOTE: we must use the > sign and not >= in case
the user enters a value of zero for q, which would
result in zero-divide.}
end;
end;
Bvec[i] := s;
end;
writeln(‘Least squares solution’);
for j := 1 to nCo1 do
begin
write(Bvec[j]:12,‘’);
if j = 5 * (j div 5) then writeln;
end;
writeln;
s := resids(nRow, nCo1, A, Y, Bvec, true);
end; {if q>=0.0}
until q<0.0; {this is how we exit from the procedure}
end {alg02.pas == svdlss};

In the above code the residual sum of squares is computed in the separate procedure resids.pas.
In alg02.pas, I have not included step numbers because the present code is quite different from the
original algorithm.