svd.sa
Generated by gen_html_sa_files from ICSI. Contact gomes@icsi.berkeley.edu for details
---------------------------> Sather 1.1 source file <--------------------------
-- COPYRIGHT NOTICE: This code is provided WITHOUT ANY WARRANTY
-- and is subject to the terms of the SATHER LIBRARY GENERAL PUBLIC
-- LICENSE contained in the file: Sather/Doc/License of the
-- Sather distribution. The license is also available from ICSI,
-- 1947 Center St., Suite 600, Berkeley CA 94704, USA.
-- This class had calls to a c version that seems to no longer exist!
-- I commented that stuff out and added the test to test-oblig.sa
-- (gomes Nov21/95)
-- *************************************************************
-- Converted to use the new matrix and vector classes July 96
-- The test class for SVD seems to generate results that are slightly
-- off from the expected results in 0.2.
-- Bo Hagstedt<bo.hagstedt@mailbox.swipnet.se> looked at the code
-- and found that it was good (and fixed several errors).
-- The IF statments protect accesses outside the matrices when
-- the loop index i points to the last row column as the sub loop will never
-- start when the start index is greater han the stop indexand a sub loop uses
-- i+1 as its index.
-- Thanks also to Alex Cozzi who looked at this code earlier.
-- gomes@icsi.berkeley.edu
class NR_SVD
class NR_SVD is
-- Singular value decomposition derived from Numerical Recipes.
-- This algorithm is known to have a problem with an obscure case (but
-- works almost all of the time). Better algorithms are being worked on.
pythag(a,b:FLT):FLT is
-- Square root of a^2+b^2 without destructive overflow or underflow.
res:FLT;
at:FLT:=a.abs; bt:FLT:=b.abs; ct:FLT;
if at>bt then ct:=bt/at; res:=at*(1.0+ct*ct).sqrt;
elsif bt/=0.0 then ct:=at/bt; res:=bt*(1.0+ct*ct).sqrt;
end; -- if
return res;
end; -- pythag
svd(a:MAT, w:VEC, v:MAT)
pre a.nr>=a.nc and w.dim=a.nc and v.nr=a.nc and v.nc=a.nc is
-- Computes the singular value decomposition of a. When finished,
-- a w v^T equals the old a. Must have a.nr>=a.nc.
-- Based on Numerical Recipes in C, p. 68.
m:INT:=a.nr; n:INT:=a.nc;
-- Householder reduction to bidiagonal form (should be separate).
rv1:VEC:=rv1.create(n);
g,scale,anorm,c,f,h,s,x,y,z:FLT;
i,its,j,jj,k,l,nm:INT; flag:BOOL;
-- Householder reduction to bidiagonal form
i:=0; loop until!(i>=n);
l:=i+1; rv1[i]:=scale*g; g:=0.0; s:=0.0; scale:=0.0;
if i<m then
k:=i; loop until!(k>=m); scale:=scale+a[k,i].abs; k:=k+1 end;
if scale/=0.0 then
k:=i; loop until!(k>=m);
a[k,i]:=a[k,i]/scale; s:=s+a[k,i]*a[k,i]; k:=k+1
end; -- loop
f:=a[i,i];
if f<0.0 then g:=s.sqrt else g:=-s.sqrt end;
h:=f*g-s; a[i,i]:=f-g;
j:=l; loop until!(j>=n);
s:=0.0;
k:=i; loop until!(k>=m); s:=s+a[k,i]*a[k,j]; k:=k+1 end;
f:=s/h;
k:=i; loop until!(k>=m); a[k,j]:=a[k,j]+f*a[k,i];k:=k+1 end;
j:=j+1
end; -- loop
k:=i; loop until!(k>=m); a[k,i]:=a[k,i]*scale; k:=k+1 end;
end; -- if
end; -- if
w[i]:=scale*g; g:=0.0; s:=0.0; scale:=0.0;
if i<m and i/=n-1 then
k:=l; loop until!(k>=n); scale:=scale+a[i,k].abs; k:=k+1 end;
if scale/=0.0 then
k:=l; loop until!(k>=n);
a[i,k]:=a[i,k]/scale; s:=s+a[i,k]*a[i,k]; k:=k+1
end; -- loop
f:=a[i,l];
if f<0.0 then g:=s.sqrt else g:=-s.sqrt end;
h:=f*g-s; a[i,l]:=f-g;
k:=l; loop until!(k>=n);
rv1[k]:=a[i,k]/h;
k:=k+1
end;
j:=l; loop until!(j>=m);
s:=0.0;
k:=l; loop until!(k>=n); s:=s+a[j,k]*a[i,k]; k:=k+1 end;
k:=l; loop until!(k>=n); a[j,k]:=a[j,k]+s*rv1[k];k:=k+1 end;
j:=j+1
end; -- loop
k:=l; loop until!(k>=n); a[i,k]:=a[i,k]*scale; k:=k+1 end;
end; -- if
end; -- if
anorm:=anorm.max(w[i].abs+rv1[i].abs);
i:=i+1;
end; -- loop
-- Accumulation of right-hand transformations.
i:=n-1; loop until!(i<0);
if i<n-1 then
if g/=0.0 then
j:=l; loop until!(j>=n);
-- double division to avoid possible underflow
v[j,i]:=(a[i,j]/a[i,l])/g;
j:=j+1
end;
j:=l; loop until!(j>=n);
s:=0.0;
k:=l; loop until!(k>=n); s:=s+a[i,k]*v[k,j]; k:=k+1 end;
k:=l; loop until!(k>=n); v[k,j]:=v[k,j]+s*v[k,i]; k:=k+1 end;
j:=j+1
end; -- loop
end; -- if
j:=l; loop until!(j>=n); v[i,j]:=0.0; v[j,i]:=0.0; j:=j+1 end;
end; -- if
v[i,i]:=1.0; g:=rv1[i]; l:=i;
i:=i-1;
end; -- loop
-- Accumulation of left-hand transformations.
-- the new code is different : i=IMIN(m,n)
-- i:=n-1; loop until!(i<0);
i:=n.min(m)-1; loop until!(i<0);
l:=i+1; g:=w[i];
-- Redundant if removed by BH
j:=l; loop until!(j>=n); a[i,j]:=0.0; j:=j+1 end;
if g/=0.0 then
g:=1.0/g;
-- if i/=n-1 then BH
j:=l; loop until!(j>=n);
s:=0.0;
k:=l; loop until!(k>=m); s:=s+a[k,i]*a[k,j]; k:=k+1 end;
f:=(s/a[i,i])*g;
k:=i; loop until!(k>=m); a[k,j]:=a[k,j]+f*a[k,i]; k:=k+1 end;
j:=j+1
end; -- loop
-- end; -- if
j:=i; loop until!(j>=m); a[j,i]:=a[j,i]*g; j:=j+1 end;
else
j:=i; loop until!(j>=m); a[j,i]:=0.0; j:=j+1 end;
end; -- if
a[i,i]:=a[i,i]+1.0;
i:=i-1;
end; -- loop
-- Diagonalization of the bidiagonal form.
k:=n-1; loop until!(k<0); -- loop over singular values
--its:=0; loop until!(its>=30); -- loop over allowed iterations
-- Should read as follows in order to get the "its" error message. /BH
its:=1; loop until!(its>30); -- loop over allowed iterations
flag:=true;
l:=k; loop until!(l<0); -- test for splitting
nm:=l-1; -- rv[1] is always zero
if rv1[l].abs.flt+anorm=anorm then flag:=false; break! end;
if w[nm].abs.flt+anorm=anorm then break! end;
l:=l-1
end; -- loop
if flag then
c:=0.0; s:=1.0;
i:=l; loop until!(i>k);
f:=s*rv1[i];
-- Corrected the typo, rv1[1] => rv1[i]. /BH
rv1[i]:=c*rv1[i];
if f.abs.flt+anorm/=anorm then
g:=w[i]; h:=pythag(f,g); w[i]:=h; h:=1.0/h;
c:=g*h; s:=-f*h;
j:=0; loop until!(j>=m);
y:=a[j,nm]; z:=a[j,i];
a[j,nm]:=y*c+z*s; a[j,i]:=z*c-y*s;
j:=j+1
end; -- loop
end; -- if
i:=i+1
end; -- loop
end; -- if
z:=w[k];
if l=k then -- convergence
if z<0.0 then -- singular value is made non-negative
w[k]:=-z;
j:=0; loop until!(j>=n); v[j,k]:=-v[j,k]; j:=j+1 end;
end; -- if
break!;
end; -- if
if its=30 then #ERR+"SVD: 30 iterations, no convergence\n"; end;
-- shift from bottom 2 by 2 minor
x:=w[l]; nm:=k-1; y:=w[nm]; g:=rv1[nm]; h:=rv1[k];
f:=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y); g:=pythag(f,1.0);
if f>=0.0 then f:=((x-z)*(x+z)+h*((y/(f+g.abs))-h))/x
else f:=((x-z)*(x+z)+h*((y/(f-g.abs))-h))/x end;
-- Next QR transformation
c:=1.0; s:=1.0;
j:=l; loop until!(j>=nm+1);
i:=j+1; g:=rv1[i]; y:=w[i]; h:=s*g; g:=c*g; z:=pythag(f,h);
rv1[j]:=z; c:=f/z; s:=h/z; f:=x*c+g*s; g:=g*c-x*s;
h:=y*s; y:=y*c;
jj:=0; loop until!(jj>=n);
x:=v[jj,j]; z:=v[jj,i]; v[jj,j]:=x*c+z*s; v[jj,i]:=z*c-x*s;
jj:=jj+1
end; -- loop
z:=pythag(f,h); w[j]:=z;
if z/=0.0 then z:=1.0/z; c:=f*z; s:=h*z end;
f:=(c*g)+(s*y); x:=(c*y)-(s*g);
jj:=0; loop until!(jj>=m);
y:=a[jj,j]; z:=a[jj,i]; a[jj,j]:=y*c+z*s; a[jj,i]:=z*c-y*s;
jj:=jj+1;
end; -- loop
j:=j+1
end; -- loop
rv1[l]:=0.0; rv1[k]:=f; w[k]:=x;
its:=its+1;
end; -- loop
k:=k-1
end; -- loop
end; -- svd
end; -- class NR_SVD
class TEST_SVD
class TEST_SVD is
include TEST;
main is
class_name("SVD");
a:MAT:=MAT::create(3,4);
a.inplace_uniform_random;
svd_matrix_test(1, a);
a:=MAT::create(3,4);
a.inplace_uniform_random;
svd_matrix_test(2, a);
a:=MAT::create(4,3); a.inplace_uniform_random;
svd_matrix_test(3, a);
a:=MAT::create(5, 5); a.inplace_uniform_random;
svd_matrix_test(4, a);
a:=MAT::create(||0.0,1.0,0.0|,|0.0,1.0,1.0|,|0.0,0.0,0.0||);
svd_matrix_test(5, a); -- the bad example for Numerical Recipes
svd_back_sub_test(1,
MAT::create(||1.0,0.0,0.0|,|0.0,2.0,0.0|,|0.0,0.0,3.0||),
VEC::create(|1.0,0.0,0.0|));
svd_back_sub_test(2,
MAT::create(||1.0,0.0,0.0|,|0.0,2.0,0.0|,|0.0,0.0,3.0||),
VEC::create(|1.0,1.0,1.0|));
svd_back_sub_test(3,
MAT::create(||1.0,2.0,3.0|,|2.0,2.0,3.0|,|3.0,2.0,3.0||),
VEC::create(|1.0,1.0,1.0|));
svd_back_sub_test(4,
MAT::create(||1.0,2.0,3.0|,|2.0,2.0,3.0|,|3.0,2.0,3.0||),
VEC::create(|1.0,2.0,3.0|));
svd_back_sub_test(5,
MAT::create(||1.0,2.0,3.0,4.0,5.0|,|2.0,3.0,4.0,5.0,1.0|
,|3.0,4.0,5.0,1.0,2.0|,|4.0,5.0,1.0,2.0,3.0|,|5.0,1.0,2.0,3.0,4.0||),
VEC::create(|1.0,1.0,1.0,1.0,1.0|));
svd_back_sub_test(6,
MAT::create(||1.0,2.0,3.0,4.0,5.0|,|2.0,3.0,4.0,5.0,1.0|
,|3.0,4.0,5.0,1.0,2.0|,|4.0,5.0,1.0,2.0,3.0|,|5.0,1.0,2.0,3.0,4.0||),
VEC::create(|1.0,2.0,3.0,4.0,5.0|));
svd_back_sub_test(7,
MAT::create(||1.4, 2.1, 2.1, 7.4, 9.6|
, |1.6, 1.5, 1.1, 0.7, 5.0|, |3.8, 8.0, 9.6, 5.4, 8.8|
, |4.6, 8.2, 8.4, 0.4, 8.0|, |2.6, 2.9, 0.1, 9.6, 7.7||),
VEC::create(|1.1, 1.6, 4.7, 9.1, 0.1|));
lfm:MAT:=MAT::create(||1.0,2.0|,|3.0,4.0|,|5.0,6.0||);
fl_lfi:FLIST{VEC}:=FLIST{VEC}::create;
fl_lfi:=fl_lfi.push(VEC::create(|1.0,0.0|));
fl_lfi:=fl_lfi.push(VEC::create(|1.0,1.0|));
fl_lfi:=fl_lfi.push(VEC::create(|0.0,1.0|));
fl_lfi:=fl_lfi.push(VEC::create(|2.0,1.0|));
lfi:ARRAY{VEC} := fl_lfi.array;
fl_lfo:FLIST{VEC}:=FLIST{VEC}::create;
loop
fl_lfo:=fl_lfo.push(lfm.times_vec(lfi.elt!));
end;
lfo: ARRAY{VEC} := fl_lfo.array;
loop
#OUT +"lfo: "+lfo.elt!.str+"\n";
end;
--d test("to_linear_fit_of",
--d lfm.copy.inplace_zero.to_linear_fit_of(lfi,lfo).str, lfm.str);
fl_lfoa:FLIST{VEC}:=FLIST{VEC}::create;
loop
fl_lfoa:=fl_lfoa.push(lfo.elt!+(VEC::create(|1.0,2.0,3.0|)));
end;
lfoa:ARRAY{VEC} := fl_lfoa.array;
lfma:MAT:=MAT::create
(||1.0,2.0,1.0|,|3.0,4.0,2.0|,|5.0,6.0,3.0||);
--d test("to_affine_fit_of",
--d lfma.copy.inplace_zero.inplace_affine_fit_of(lfi,lfoa).str, lfma.str);
--d wt:FLIST{FLT}:=FLIST{FLT}::create;
--d wt:=wt.push(1.0); wt:=wt.push(1.0); wt:=wt.push(1.0); wt:=wt.push(1.0);
--d lfmcopy ::= lfm.copy;
--d lfmcopy.inplace_zero;
--d lfmcopy.inplace_weighted_linear_fit_of(lfi,lfo,wt);
--d test("to_weighted_linear_fit_of",lfmcopy.str,lfm.str);
--d test("to_weighted_affine_fit_of",
--d lfma.copy.inplace_zero.inplace_weighted_affine_fit_of(lfi,lfoa,wt).str
--d , lfma.str);
finish;
end;
svd_matrix_test(n:INT, a:MAT) is
-- Test the singular value decomposition of `a' labelling the test
-- by `n'.
u:MAT:=u.create(a.nr.max(a.nc),a.nc);
v:MAT:=v.create(a.nc, a.nc); w:VEC:=w.create(a.nc);
a.svd_in(u,w,v);
u.inplace_times_diagonal(w);
v.inplace_trans;
tmp:MAT:=u.copy;
tmp.inplace_arg_times_arg(u,v);
acpy::= a.copy;
acpy.inplace_zero;
acpy.inplace_portion_of_arg(tmp);
cmp:MAT:=acpy;
unchecked_test("svd_in "+n, cmp.str, a.str);
end; -- svd_test_matrix
svd_back_sub_test(n:INT, a:MAT, b:VEC) is
-- Test `svd_back_sub' on the problem `a.x=b', label with `n'.
size:INT:=a.nr;
u:MAT:=u.create(size,size); v:MAT:=v.create(size,size);
w:VEC:=w.create(size); a.svd_in(u,w,v);
wmax:FLT:=w.max_value; wmin:FLT:=wmax*(0.000001);
k:INT:=0; loop until!(k=size);
if w[k]<=wmin then w[k]:=0.0 end; k:=k+1
end; -- loop
x:VEC:=VEC::create(size);
a.svd_back_sub(u,w,v,b,x);
c:VEC:=a.times_vec(x);
unchecked_test("svd_back_sub "+n, c.str, b.str);
end; -- svd_back_sub_test
end; -- class TEST_SVD