One-sided Jacobi implementation of SVD

Question

I'm trying to write a simple implementation of Singular Value Decomposition (SVD). I'm using the one-sided Jacobi algorithm since it seemed like a simple one. The algorithm is described here and there's a simple Matlab code for that here (exercise 4). I've implemented the same code and it works ok and by that I mean that SVD(A) = U * S * V' (or any other notation that's used) and for some matrices the result is the same as what Matlab's SVD produces (except this function doesn't sort the singular values). But my problem is that when the matrix A has a 0 singular value, U and V are no longer unitary matrices as they should be!

Is there a way to update this algorithm so that it works for cases with 0 singular value too? If not, is there another algorithm for SVD that's as easy to implement? Any help is appreciated.

Here's my Matlab code that's basically same as the one in the link above, just completed and with minor changes.

function [U,S,V] = jacobi_SVD(A)
  TOL=1.e-4;
  n=size(A,1);
  U=A';
  V=eye(n);
  converge=TOL+1;
  while converge>TOL
  converge=0;
  for i=1:n-1
    for j=i+1:n
      % compute [alpha gamma;gamma beta]=(i,j) submatrix of U*U'
      alpha=sumsqr(U(i, :)); 
      beta=sumsqr(U(j, :));
      gamma=sum(U(i, :).* U(j, :));
      converge=max(converge,abs(gamma)/sqrt(alpha*beta));

      % compute Jacobi rotation that diagonalizes 
      %    [alpha gamma;gamma beta]
      zeta=(beta-alpha)/(2*gamma);
      t=sign(zeta)/(abs(zeta)+sqrt(1+zeta^2));
      c= 1.0 / (sqrt(1 + t * t));
      s= c * t;

      % update columns i and j of U
      t=U(i, :);
      U(i, :)=c*t-s*U(j, :);
      U(j, :)=s*t+c*U(j, :);

      % update matrix V of right singular vectors
      t=V(i, :);
      V(i, :)=c*t-s*V(j, :);
      V(j, :)=s*t+c*V(j, :);
    end
  end
end

% the singular values are the norms of the columns of U
% the left singular vectors are the normalized columns of U
for j=1:n
  singvals(j)=norm(U(j, :));
  U(j, :)=U(j, :)/singvals(j);
end
S=diag(singvals);
U = U';
V = V'; %return V, not V'
end

Unknown Unknown · Accepted Answer · 2018-06-07T16:18:24

I can't get your code to run because on when you evaluate alpha and beta you have sumsq which you haven't defined. I found some simple code available on the Matlab website. which uses the QR decomposition (Gram-schmidt).

function [u,s,v] = svdsim(a,tol)
%SVDSIM  simple SVD program
%
% A simple program that demonstrates how to use the
% QR decomposition to perform the SVD of a matrix.
% A may be rectangular and complex.
%
% usage: [U,S,V]= SVDSIM(A)
%     or      S = SVDSIM(A)
%
% with A = U*S*V' , S>=0 , U'*U = Iu  , and V'*V = Iv
%
% The idea is to use the QR decomposition on A to gradually "pull" U out from
% the left and then use QR on A transposed to "pull" V out from the right.
% This process makes A lower triangular and then upper triangular alternately.
% Eventually, A becomes both upper and lower triangular at the same time,
% (i.e. Diagonal) with the singular values on the diagonal.
%
% Matlab's own SVD routine should always be the first choice to use,
% but this routine provides a simple "algorithmic alternative"
% depending on the users' needs.
% 
%see also: SVD, EIG, QR, BIDIAG, HESS
%

% Paul Godfrey
% October 23, 2006

if ~exist('tol','var')
   tol=eps*1024;
end

%reserve space in advance
sizea=size(a);
loopmax=100*max(sizea);
loopcount=0;

% or use Bidiag(A) to initialize U, S, and V
u=eye(sizea(1));
s=a';
v=eye(sizea(2));

Err=realmax;
while Err>tol & loopcount<loopmax ;
%   log10([Err tol loopcount loopmax]); pause
    [q,s]=qr(s'); u=u*q;
    [q,s]=qr(s'); v=v*q;

% exit when we get "close"
    e=triu(s,1);
    E=norm(e(:));
    F=norm(diag(s));
    if F==0, F=1;end
    Err=E/F;
    loopcount=loopcount+1;
end
% [Err/tol loopcount/loopmax]

%fix the signs in S
ss=diag(s);
s=zeros(sizea);
for n=1:length(ss)
    ssn=ss(n);
    s(n,n)=abs(ssn);
    if ssn<0
       u(:,n)=-u(:,n);
    end
end

if nargout<=1
   u=diag(s);
end

return

Typically you don't actually have zeros in my experience. Numerical precision leaves you with something approximately near that, for instance, the following. If I want to create a 5 x 5 matrix but have it be rank 3 I can do the following.

A = randn(5,3)*randn(5,3);
[U,S,Vt] = svdsim(A,1e-8);


S =

    6.3812         0         0         0         0
         0    2.0027         0         0         0
         0         0    1.0240         0         0
         0         0         0    0.0000         0
         0         0         0         0    0.0000

Now it simply looks like you have zeros. But if you look closer.

format long 
>> S(4,4)

ans =

     3.418057860623250e-16


   S(5,5)

ans =

     9.725444388260210e-17

I will note that this is machine epsilon and for all intensive purposes it is 0.

One-sided Jacobi implementation of SVD

1 Answers