Artifact 56bf3cdf0698dccb3701783fded001eee12a33bd6fd7558c245ac6a9d2c9b917:
- Executable file
r37/packages/linalg/svd.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 16509) [annotate] [blame] [check-ins using] [more...]
- Executable file
r38/packages/linalg/svd.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 16509) [annotate] [blame] [check-ins using]
module svd; %**********************************************************************% % % % Computation of the Singular Value Decomposition of dense matrices % % containing numeric entries. Uses specific rounded number routines to % % speed things up. % % % % Author: Matt Rebbeck, June 1994. % % % % The algorithm was taken from "Linear Algebra" - J.H.Wilkinson % % & C. Reinsch % % % %**********************************************************************% symbolic smacro procedure my_minus(u); % % Efficiently performs reval({'minus,u}). % if atom u then {'minus,u} else if car u = 'minus then cadr u else {'minus,u}; symbolic procedure svd(A); % % Computation of the singular values and complete orthogonal % decomposition of a real rectangular matrix A. % % A = tp(U) diag(q) V, U tp(U) = V tp(V) = I, % % and q contains the singular values along the diagonal. % (tp => transpose). % % Algorithm taken from "Linear Algebra" % J.H.Wilkinson & C.Reinsch % begin scalar ee,U,V,g,x,eps,tolerance,q,s,f,h,y,test_f_splitting, cancellation,test_f_convergence,convergence,c,z,denom,q_mat, I_rounded_turned_on,trans_done; integer i,j,k,l,l1,m,n; trans_done := I_rounded_turned_on := nil; if not !*rounded then << on rounded; I_rounded_turned_on := t; >>; if not matrixp(A) then rederr "Error in svd: non matrix input."; % The value of eps can be decreased to increase accuracy. % As usual, doing this will slow things down (and vice versa). % It should not be made smaller than the value of rd!-tolerance!*. eps := get_num_part(my_reval({'times,1.5,{'expt,10,-8}})); tolerance := get_num_part(my_reval({'expt,10,-31})); % Algorithm requires m >= n. If this is not the case then transpose % the input and swap U and V in the output (as A = tp(U) diag(q) V % but tp(A) = tp(V) diag(q) U ). if row_dim(A) < column_dim(A) then << A := algebraic tp(A); trans_done := t; >>; m := row_dim(A); n := column_dim(A); U := rd_copy_mat(A); V := mkmatrix(n,n); ee := mkvect(n); q := mkvect(n); % Householder's reduction to bidiagonal form: g := x := 0; for i:=1:n do << putv(ee,i,g); s := 0; l := i+1; for j:=i:m do s := specrd!:plus(s,specrd!:expt(getmat(U,j,i),2)); if get_num_part(s) < tolerance then g := 0 else << f := getmat(U,i,i); if get_num_part(f)<0 then g := specrd!:sqrt(s) else g := my_minus(specrd!:sqrt(s)); h := specrd!:plus(specrd!:times(f,g),my_minus(s)); setmat(U,i,i,specrd!:plus(f,my_minus(g))); for j:=l:n do << s := 0; for k:=i:m do s := specrd!:plus(s,specrd!:times(getmat(U,k,i), getmat(U,k,j))); f := specrd!:quotient(s,h); for k:=i:m do setmat(U,k,j,specrd!:plus(getmat(U,k,j), specrd!:times(f,getmat(U,k,i)))); >>; >>; putv(q,i,g); s := 0; for j:=l:n do s := specrd!:plus(s,specrd!:expt(getmat(U,i,j),2)); if get_num_part(s) < tolerance then g := 0 else << f := getmat(U,i,i+1); if get_num_part(f)<0 then g := specrd!:sqrt(s) else g := my_minus(specrd!:sqrt(s)); h := specrd!:plus(specrd!:times(f,g),my_minus(s)); setmat(U,i,i+1,specrd!:plus(f,my_minus(g))); for j:=l:n do putv(ee,j,specrd!:quotient(getmat(U,i,j),h)); for j:=l:m do << s := 0; for k:=l:n do s := specrd!:plus(s,specrd!:times(getmat(U,j,k), getmat(U,i,k))); for k:=l:n do setmat(U,j,k,specrd!:plus(getmat(U,j,k), specrd!:times(s,getv(ee,k)))); >>; >>; y := specrd!:plus(abs(get_num_part(getv(q,i))), abs(get_num_part(getv(ee,i)))); if get_num_part(y) > get_num_part(x) then x := y; >>; % Accumulation of right hand transformations: for i:=n step -1 until 1 do << if get_num_part(g) neq 0 then << h := specrd!:times(getmat(U,i,i+1),g); for j:=l:n do setmat(V,j,i,specrd!:quotient(getmat(U,i,j),h)); for j:=l:n do << s := 0; for k:=l:n do s := specrd!:plus(s,specrd!:times(getmat(U,i,k), getmat(V,k,j))); for k:=l:n do setmat(V,k,j,specrd!:plus(getmat(V,k,j), specrd!:times(s,getmat(V,k,i)))); >>; >>; for j:=l:n do << setmat(V,i,j,0); setmat(V,j,i,0); >>; setmat(V,i,i,1); g := getv(ee,i); l := i; >>; % Accumulation of left hand transformations: for i:=n step -1 until 1 do << l := i+1; g := getv(q,i); for j:=l:n do setmat(U,i,j,0); if get_num_part(g) neq 0 then << h := specrd!:times(getmat(U,i,i),g); for j:=l:n do << s := 0; for k:=l:m do s := specrd!:plus(s,specrd!:times(getmat(U,k,i), getmat(U,k,j))); f := specrd!:quotient(s,h); for k:=i:m do setmat(U,k,j,specrd!:plus(getmat(U,k,j), specrd!:times(f,getmat(U,k,i)))); >>; for j:=i:m do setmat(U,j,i,specrd!:quotient(getmat(U,j,i),g)); >> else for j:=i:m do setmat(U,j,i,0); setmat(U,i,i,specrd!:plus(getmat(U,i,i),1)); >>; % Diagonalisation of the bidiagonal form: eps := get_num_part(specrd!:times(eps,x)); test_f_splitting := t; k := n; while k>=1 do << convergence := nil; if test_f_splitting then << l := k; test_f_convergence := cancellation := nil; while l>=1 and not (test_f_convergence or cancellation) do << if abs(get_num_part(getv(ee,l))) <= eps then test_f_convergence := t else if abs(get_num_part(getv(q,l-1))) <= eps then cancellation := t else l := l-1; >>; >>; % Cancellation of e[l] if l>1: if not test_f_convergence then << c := 0; s := 1; l1 := l-1; i := l; while i<=k and not test_f_convergence do << f := specrd!:times(s,getv(ee,i)); putv(ee,i,specrd!:times(c,getv(ee,i))); if abs(get_num_part(f)) <= eps then test_f_convergence := t else << g := getv(q,i); h := specrd!:sqrt(specrd!:plus(specrd!:times(f,f), specrd!:times(g,g))); putv(q,i,h); c := specrd!:quotient(g,h); s := specrd!:quotient(my_minus(f),h); for j:=1:m do << y := getmat(U,j,l1); z := getmat(U,j,i); setmat(U,j,l1,specrd!:plus(specrd!:times(y,c), specrd!:times(z,s))); setmat(U,j,i,specrd!:difference(specrd!:times(z,c), specrd!:times(y,s))); >>; i := i+1; >>; >>; >>; z := getv(q,k); if l = k then convergence := t; if not convergence then << % Shift from bottom 2x2 minor: x := getv(q,l); y := getv(q,k-1); g := getv(ee,k-1); h := getv(ee,k); f := specrd!:quotient(specrd!:plus(specrd!:times( specrd!:plus(y,my_minus(z)),specrd!:plus(y,z)), specrd!:times(specrd!:plus(g,my_minus(h)), specrd!:plus(g,h))),specrd!:times( specrd!:times(2,h),y)); g := specrd!:sqrt(specrd!:plus(specrd!:times(f,f),1)); % Needed to change < here to <=. if get_num_part(f)<=0 then denom := specrd!:plus(f,my_minus(g)) else denom := specrd!:plus(f,g); f := specrd!:quotient(specrd!:plus(specrd!:times( specrd!:plus(x,my_minus(z)),specrd!:plus(x,z)), specrd!:times(h,specrd!:quotient(y, specrd!:plus(denom,my_minus(h))))),x); % Next QR transformation: c := s := 1; for i:=l+1:k do << g := getv(ee,i); y := getv(q,i); h := specrd!:times(s,g); g := specrd!:times(c,g); z := specrd!:sqrt(specrd!:plus(specrd!:times(f,f), specrd!:times(h,h))); putv(ee,i-1,z); c := specrd!:quotient(f,z); s := specrd!:quotient(h,z); f := specrd!:plus(specrd!:times(x,c),specrd!:times(g,s)); g := specrd!:plus(specrd!:times(my_minus(x),s), specrd!:times(g,c)); h := specrd!:times(y,s); y := specrd!:times(y,c); for j:=1:n do << x := getmat(V,j,i-1); z := getmat(V,j,i); setmat(V,j,i-1,specrd!:plus(specrd!:times(x,c), specrd!:times(z,s))); setmat(V,j,i,specrd!:difference(specrd!:times(z,c), specrd!:times(x,s))); >>; z := specrd!:sqrt(specrd!:plus(specrd!:times(f,f), specrd!:times(h,h))); putv(q,i-1,z); c := specrd!:quotient(f,z); s := specrd!:quotient(h,z); f := specrd!:plus(specrd!:times(c,g),specrd!:times(s,y)); x := specrd!:plus(specrd!:times(my_minus(s),g), specrd!:times(c,y)); for j:=1:m do << y := getmat(U,j,i-1); z := getmat(U,j,i); setmat(U,j,i-1,specrd!:plus(specrd!:times(y,c), specrd!:times(z,s))); setmat(U,j,i,specrd!:difference(specrd!:times(z,c), specrd!:times(y,s))); >>; >>; putv(ee,l,0); putv(ee,k,f); putv(q,k,x); >> else % convergence: << if get_num_part(z)<0 then << % q[k] is made non-negative: putv(q,k,my_minus(z)); for j:=1:n do setmat(V,j,k,my_minus(getmat(V,j,k))); >>; k := k-1; >>; >>; q_mat := q_to_diag_matrix(q); if I_rounded_turned_on then off rounded; if trans_done then return {'list,algebraic tp V,q_mat,algebraic tp U} else return {'list,algebraic tp U,q_mat,algebraic tp V}; end; flag('(svd),'opfn); % To make it available from algebraic (user) mode. symbolic procedure q_to_diag_matrix(q); % % Converts q (a vector) to a diagonal matrix with the elements of % q on the diagonal. % begin scalar q_mat; integer i,sq_dim_q; sq_dim_q := upbv(q); q_mat := mkmatrix(sq_dim_q,sq_dim_q); for i:=1:sq_dim_q do setmat(q_mat,i,i,getv(q,i)); return q_mat; end; symbolic procedure pseudo_inverse(in_mat); % % Also known as the Moore-Penrose Inverse. % % Given the singular value decomposition A := tp(U) diag(q) V % the pseudo inverse A^(-1) is defined as % % A^(-1) = tp(V) (diag(q))^(-1) U. % % NB: this can be quite handy as we can take the inverse of non % square matrices (A * pseudo_inverse(A) = identity). % begin scalar psu_inv,svd_list; svd_list := svd(in_mat); psu_inv := algebraic (tp(third svd_list)*(1/second svd_list)*first svd_list); return psu_inv; end; flag('(pseudo_inverse),'opfn); symbolic procedure rd_copy_mat(A); % % Creates a copy of the input matrix and returns it aswell as % reval-ing each elt to get them in !:rd!: form; % begin scalar C; integer row_dim,column_dim; row_dim := first size_of_matrix(A); column_dim := second size_of_matrix(A); C := mkmatrix(row_dim,column_dim); for i:=1:row_dim do << for j:=1:column_dim do << setmat(C,i,j,my_reval(getmat(A,i,j))); >>; >>; return C; end; % % All computation is done with rounded mode on and with all numbers % in !:rd!: form. The following specrd!: functions makes the algebraic % computation of these numbers very efficient. % symbolic procedure specrd!:times(u,v); begin scalar negsign; u := add_minus(u); v := add_minus(v); if eqcar(u,'minus) then << u := cadr u; negsign := t>>; if eqcar(v,'minus) then << v := cadr v; negsign := not negsign>>; if atom u then u := mkround float u; if atom v then v := mkround float v; return if negsign then list('minus,rd!:times(u,v)) else rd!:times(u,v); end; symbolic procedure specrd!:quotient(u,v); begin scalar negsign; u := add_minus(u); v := add_minus(v); if eqcar(u,'minus) then << u := cadr u; negsign := t>>; if eqcar(v,'minus) then << v := cadr v; negsign := not negsign>>; if atom u then u := mkround float u; if atom v then v := mkround float v; return if negsign then list('minus,rd!:quotient(u,v)) else rd!:quotient(u,v); end; symbolic procedure specrd!:expt(u,v); begin if (u = '(!:rd!: . 0.0) or u = 0) then return '(!:rd!: . 0.0); if eqcar(u,'minus) then u := ('!:rd!: . -cdadr u); if eqcar(v,'minus) then v := ('!:rd!: . -cdadr v); if atom u then u := mkround float u; if atom v then v := mkround float v; return rdexpt!*(u,v); end; symbolic procedure specrd!:sqrt(u); specrd!:expt(u,0.5); symbolic procedure specrd!:plus(u,v); begin scalar negsign; negsign := 0; u := add_minus(u); v := add_minus(v); if eqcar(u,'minus) then << u := cadr u; negsign := 1>>; if eqcar(v,'minus) then << v := cadr v; negsign := negsign +2>>; if atom u then u := mkround float u; if atom v then v := mkround float v; return if negsign = 0 then rd!:plus(u,v) else if negsign = 3 then list('minus,rd!:plus(u,v)) else if negsign =2 then rd!:difference (u,v) else rd!:difference(v,u); end; symbolic procedure specrd!:difference(u,v); begin scalar negsign; negsign := 0; u := add_minus(u); v := add_minus(v); if eqcar(u,'minus) then << u := cadr u; negsign := 1>>; if eqcar(v,'minus) then << v := cadr v; negsign := negsign +2>>; if atom u then u := mkround float u; if atom v then v := mkround float v; return if negsign = 0 then rd!:difference(u,v) else if negsign = 3 then list('minus,rd!:difference(u,v)) else if negsign =2 then rd!:plus (u,v) else list('minus,rd!:plus(v,u)); end; symbolic procedure add_minus(u); % % Things like (!:rd!: . -0.12345) can cause problems as negsign does % not notice the negative. This function converts that to % {'minus,(!:rd!: . 0.12345)}. Unfortunately it slows things down but % it works. % begin if atom u then return u else if car u = '!:rd!: and cdr u >= 0 then return u else if car u = '!:rd!: and cdr u < 0 then return {'minus,('!:rd!: . abs(cdr u))} else if car u = 'minus and numberp cadr u then return u else if car u = 'minus and cdadr u < 0 then return ('!:rd!: . abs(cdadr u)) else if car u = 'minus then return u else if cdr u < 0 then return {'minus,('!:rd!: . abs(cdr u))} else return u; end; endmodule; % svd. end;