奇异值分解的单边Jacobi实现

Question

我正在尝试编写一个简单的奇异值分解(SVD)实现。我使用单边Jacobi算法，因为它看起来很简单。该算法在here中进行了描述，该here有一个简单的Matlab代码(练习4)。我已经实现了相同的代码，它工作得很好，我的意思是SVD(A) =U*S* V‘(或使用的任何其他表示法)，对于某些矩阵，结果与Matlab的SVD产生的结果相同(除了这个函数不对奇异值进行排序)。但是我的问题是，当矩阵A有一个0的奇异值时，U和V不再是它们应该是的么正矩阵！

有没有办法更新这个算法，让它也适用于奇异值为0的情况？如果没有，有没有另一种同样容易实现的SVD算法？任何帮助都是非常感谢的。

这是我的Matlab代码，它与上面链接中的代码基本相同，只是完成了一些微小的更改。

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

Answer 1

我不能让你的代码运行，因为当你评估alpha和beta时，你得到了你还没有定义的sumsq。我在Matlab website.上找到了一些使用QR分解(格拉姆-施密特)的简单代码。

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

在我的经验中，通常情况下你实际上没有零。数值精度给你留下了大致接近的东西，例如，下面的内容。如果我想创建一个5 x 5的矩阵，但它的秩是3，我可以这样做。

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

现在它看起来就像是你有零。但如果你仔细观察。

format long 
>> S(4,4)

ans =

     3.418057860623250e-16


   S(5,5)

ans =

     9.725444388260210e-17

我要注意的是，这是机器epsilon，对于所有密集的目的，它是0。