Разделяй и властвуй SVD в MATLAB

Я пытаюсь реализовать разделяй и властвуй 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

matlab numeric numerical-methods numerical-computing
person dxdydz    schedule 31.03.2020    source источник
comment
Должно ли это B1 = B(1:k-1,1:k); быть B1 = B(1:k,1:k);? Кроме того, ваш код неверен, например. для B = [-1] (матрица S должна содержать только положительные значения, вам нужно сделать либо U, либо V знаком B, а S - абсолютным значением B для матриц 1x1)   -  person chtz    schedule 31.03.2020
comment
какой размер B вы используете? Я не смог воспроизвести вашу ошибку   -  person Yuval Harpaz    schedule 05.04.2020
comment
@Yuval Harpaz, я использовал матрицу 6 на 6   -  person dxdydz    schedule 05.04.2020
comment
теперь я получаю сообщение об ошибке, а как насчет комментария @chtz о B1 = B (1: k, 1: k)? если вы удалите все -1 из своего кода, вы не получите ошибки. И если вы этого не сделаете, вы никогда не используете B (k, k)   -  person Yuval Harpaz    schedule 06.04.2020

Ответы (1)


arrow_upward
0
arrow_downward

Имея ограниченные математические знания, вам нужно немного помочь мне, так как я не могу судить, правилен ли подход с математической точки зрения (без какой-либо теории;)). В любом случае, я даже не смог воспроизвести ошибку, например, с этой матрицей, которую 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
person max    schedule 05.04.2020