Разделяй и властвуй СВД в MATLAB - PullRequest
Новая
       8

Разделяй и властвуй СВД в MATLAB

1 голос
/ 31 марта 2020

Я пытаюсь реализовать разделение и завоевание SVD верхней двухдиагональной матрицы B, но мой код не работает. Ошибка:

"Невозможно выполнить назначение, потому что размер левой стороны 3 на 3, а размер правой стороны 2 на 2.
V_bar (1: k , 1: k) = V1; "

Может кто-нибудь помочь мне это исправить? Спасибо. 

function [U,S,V] = DivideConquer_SVD(B)
[m,n] = size(B);
k = floor(m/2);
if k == 0
   U = 1;
   V = 1;
   S = B;
   return;
else
    % Divide the input matrix
    alpha = B(k,k);
    beta = B(k,k+1);
    e1 = zeros(m,1);
    e2 = zeros(m,1);
    e1(k) = 1;
    e2(k+1) = 1;
    B1 = B(1:k-1,1:k);
    B2 = B(k+1:m,k+1:m);
    %recursive computations
    [U1,S1,V1] = DivideConquer_SVD(B1);
    [U2,S2,V2] = DivideConquer_SVD(B2);
    U_bar = zeros(m);
    U_bar(1:k-1,1:k-1) = U1;
    U_bar(k,k) = 1;
    U_bar((k+1):m,(k+1):m) = U2;
    D = zeros(m);
    D(1:k-1,1:k) = S1;
    D((k+1):m,(k+1):m) = S2;
    V_bar = zeros(m);
    V_bar(1:k,1:k) = V1;
    V_bar((k+1):m,(k+1):m) = V2;
    u = alpha*e1'*V_bar + beta*e2'*V_bar;
    u = u';
    D_tilde = D*D + u*u';
    % compute eigenvalues and eigenvectors of D^2+uu'
    [L1,Q1] = eig(D_tilde);
    eigs = diag(L1);
    S = zeros(m,n)
    S(1:(m+1):end) = eigs
    U_tilde = Q1;
    V_tilde = Q1;
    %Compute eigenvectors of the original input matrix T
    U = U_bar*U_tilde;
    V = V_bar*V_tilde;
    return;
end

1 Ответ

0 голосов
/ 05 апреля 2020

Имея ограниченные математические знания, вы должны помочь мне немного больше - поскольку я не могу судить, правильн ли подход математически (без теории;)). В любом случае, я даже не смог воспроизвести ошибку, например, с помощью этой матрицы, которую MathWorks использует для иллюстрации их LU-факторизации матрицы 

A = [10 -7 0
     -3  2 6
      5 -1 5];

Поэтому я попытался немного структурировать ваш код и дал несколько подсказок. Расширьте это, чтобы сделать ваш код более понятным для тех людей (как я), которые не слишком знакомы с разложением матриц. 

function [U,S,V] = DivideConquer_SVD(B)
% m x n matrix
[m,n] = size(B);

k = floor(m/2);
if k == 0
    disp('if') % for debugging
   U = 1;
   V = 1;
   S = B;
%    return; % net necessary as you don't do anything afterwards anyway
else
    disp('else') % for debugging
    % Divide the input matrix
    alpha = B(k,k); % element on diagonal
    beta = B(k,k+1); % element on off-diagonal
    e1 = zeros(m,1);
    e2 = zeros(m,1);
    e1(k) = 1;
    e2(k+1) = 1;
    % divide matrix 
    B1 = B(1:k-1,1:k); % upper left quadrant
    B2 = B(k+1:m,k+1:m); % lower right quadrant
    % recusrsive function call
    [U1,S1,V1] = DivideConquer_SVD(B1);
    [U2,S2,V2] = DivideConquer_SVD(B2);

    U_bar = zeros(m);
    U_bar(1:k-1,1:k-1) = U1;
    U_bar(k,k) = 1;
    U_bar((k+1):m,(k+1):m) = U2;

    D = zeros(m);
    D(1:k-1,1:k) = S1;
    D((k+1):m,(k+1):m) = S2;

    V_bar = zeros(m);
    V_bar(1:k,1:k) = V1;
    V_bar((k+1):m,(k+1):m) = V2;

    u = (alpha*e1.'*V_bar + beta*e2.'*V_bar).'; % (little show-off tip: ' 
    % is the complex transpose operator; .' is the "normal" transpose 
    % operator. It's good practice to distinguish between them but there 
    % is no difference for real matrices anyway)


    D_tilde = D*D + u*u.';
    % compute eigenvalues and eigenvectors of D^2+uu'
    [L1,Q1] = eig(D_tilde);
    eigs = diag(L1);

    S = zeros(m,n);
    S(1:(m+1):end) = eigs;
    U_tilde = Q1;
    V_tilde = Q1;
    % Compute eigenvectors of the original input matrix T
    U = U_bar*U_tilde;
    V = V_bar*V_tilde;
%    return; % net necessary as you don't do anything afterwards anyway
end % for
end % function
...