时间戳

首页 > 解决方案 > SVD 的单面 Jacobi 实现

matlab - SVD 的单面 Jacobi 实现

问题描述

我正在尝试编写奇异值分解（SVD）的简单实现。我正在使用单边 Jacobi 算法，因为它看起来很简单。该算法在此处进行了描述，并且此处有一个简单的 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

标签： matlabsvd

解决方案

我无法让您的代码运行，因为当您评估 alpha 和 beta 时，您有未定义的 sumsq。我在 Matlab网站上找到了一些简单的代码。它使用 QR 分解 (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

通常，根据我的经验，您实际上并没有零。数值精度使您得到的结果大致接近，例如，以下。如果我想创建一个 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。

推荐阅读