chenmh_svd.m

Frequently used algorithms for numerical analysis and signal processing

function [q,U,V] = Chenmh_SVD(A)
% SVDGR  is a MATLAB translation of the Golub-Reinsch ALGOL60 procedure to
%        compute the singular value decomposition of a matrix. 
%        [q,U,V] = svdgr(A) computes two orthogonal matrices U and V
%        such that  A = U diag(q) V' 
%        Walter Gander, 14.3.2005

if nargout == 3, 
   withu = true; withv = true;
elseif nargout == 2,
    withu=true; withv = false;
else
  withu = false; withv = false;
end

tol = realmin/eps; 
[m,n] = size(A); 
U = A;

% Hoseholder's reduction to bidiagonal form
g=0; x=0; 
for i=1:n
   e(i) = g;
   
   s = norm(U(i:m,i))^2; 
   if s<tol; g=0;
   else 
     f = U(i,i); 
     if f<0; g = sqrt(s); else g = -sqrt(s); end
     h = f*g-s; U(i,i) = f-g;
     for j=(i+1):n
       s = U(i:m,i)'*U(i:m,j);
       f = s/h;
       U(i:m,j) = U(i:m,j) + f*U(i:m,i); 
     end
   end
   q(i) = g;
   
   s = norm(U(i,(i+1):n))^2;
   if s<tol; g=0;
   else
     f=U(i,i+1); 
     if f<0; g = sqrt(s); else g = -sqrt(s); end
     h = f*g-s; U(i,i+1) = f-g;
     e((i+1):n) = U(i,(i+1):n)/h;
     for j=(i+1):m
       s = U(j,(i+1):n)*U(i,(i+1):n)';
       U(j,(i+1):n) = U(j,(i+1):n) +s*e((i+1):n);
     end
   end
   y = abs(q(i))+ abs(e(i)); if y>x; x=y; end   
end %for i

% Accumulation of  right-hand transformations;
if withv; 
  for i= n:-1:1
    if g ~= 0;
      h=U(i,i+1)*g;
      V((i+1):n,i) = U(i,(i+1):n)'/h;
      for j=(i+1):n
        s = U(i,(i+1):n)*V((i+1):n,j);
        V((i+1):n,j) = V((i+1):n,j) + s*V((i+1):n,i);
      end
    end
    V(i,(i+1):n) = 0; V((i+1):n,i)=0;
    V(i,i) = 1; g = e(i); 
  end
end 


% Accumulation of  left-hand transformations;
if withu;
  for i= n:-1:1
    g = q(i);
    for j = (i+1):n, U(i,j)=0; end
    if g ~= 0;
      h = U(i,i)*g;
      for j=(i+1):n
         s = U((i+1):m,i)'*U((i+1):m,j);
         f = s/h;
         U(i:m,j) = U(i:m,j) + f*U(i:m,i);
      end
      U(i:m,i) =  U(i:m,i)/g;
    else
     U(i:m,i) = 0;
    end
    U(i,i) = U(i,i) +1;
  end % for i
end %withu

% Diagonalization of the bidiagonal form 

eps2 = eps*x;
for k= n:-1:1
  % test f splitting
  convergence = false;
  while ~convergence 
     l = k; 						     
     cancellation = false;					     
     testfconvergence = false;				     
     while (l>= 1) &  (~testfconvergence) & (~cancellation),
       if abs(e(l)) <= eps2; testfconvergence = true;	    
       elseif abs(q(l-1)) <= eps2,
              cancellation = true; 
       else l = l-1;
       end 						     
     end 							     

    % cancellation of e(l) if l>1				     
      if cancellation,		
        c=0; s=1; l1 = l-1;				     
        i= l; testf = false;				     
        while (i<=k) & (~testf) 				     
           f = s*e(i); e(i) = c*e(i); 			     
           testf = abs(f) <= eps2;			     
           if ~testf;					     
             g = q(i); h= sqrt(f*f+g*g); q(i)=h; c=g/h; s=-f/h; 
             if withu;					     
               y = U(:,l1); z = U(:,i);			     
               U(:,l1) = y*c+z*s; U(:,i) = -y*s+z*c;	     
             end % withu					     
           end % testf					     
           i=i+1;						     
        end % while i					     
      end % cancellation           
      z = q(k); convergence = l==k;
      if ~convergence
        % shift from bottom 2x2 minor
        x = q(l); y=q(k-1); g = e(k-1); h=e(k);
        f = ((y-z)*(y+z) + (g-h)*(g+h))/(2*h*y); g= sqrt(f*f+1);
        if f<0; fg = f-g; else fg = f+g; end
        f = ((x-z)*(x+z) + h*(y/fg)-h)/x;
      % next QR transformation
        c = 1; s = 1;
        for i=l+1:k
          g = e(i); y = q(i); h = s*g; g = c*g;
          z = sqrt(f*f+h*h); e(i-1) =z; c = f/z; s=h/z;
          f = x*c + g*s; g = -x*s+g*c; h = y*s; y = y*c;
          if withv
             x = V(:,i-1); z = V(:,i);
             V(:,i-1) = x*c+z*s; V(:,i) = -x*s+z*c;
          end % if withv
          z = sqrt(f*f+h*h); q(i-1) = z; c = f/z; s = h/z;
          f = c*g+s*y; x = -s*g+c*y;
          if withu, 
             y = U(:,i-1); z = U(:,i);
             U(:,i-1) = y*c+z*s; U(:,i) = -y*s+z*c;
          end %if withu
        end %for i
        e(l)=0; e(k) = f; q(k) = x;
      end % if not convergence 
    end % while ~convergence
    % convergence: 
    if z<0,
      % y(k) is made non-negative
      q(k) = -z;
      if withv, V(:,k) = -V(:,k); end
    end
end % for k


% sort singular values for compatibility with Matlab SVD 
% ascending order
[q,I] = sort(q); 
if withv, V = V(:,I); end
if withu; U = U(:,I); end

% decending order
q = q(n:-1:1)'; 
if withv, V = V(:,n:-1:1); end 
if withu, U = U(:,n:-1:1); end