Односторонняя реализация SVD Якоби

Question

Я пытаюсь написать простую реализацию разложения по сингулярным числам (SVD).Я использую односторонний алгоритм Якоби, так как он показался мне простым.Алгоритм описан здесь , и для этого есть простой код Matlab здесь (упражнение 4).Я реализовал тот же код, и он работает нормально, и я имею в виду, что SVD (A) = U * S * V '(или любая другая используемая запись), и для некоторых матриц результат такой же, как и для SVD Matlab(кроме этой функции не сортирует особые значения).Но моя проблема в том, что когда матрица 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

Я не могу запустить ваш код, потому что при оценке альфа и бета у вас есть 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

Я отмечу, что это машинный эпсилон и для всех интенсивных целей он равен 0.