MATLAB 中的分而治之 SVD答案

【问题标题】：Divide and Conquer SVD in MATLABMATLAB 中的分而治之 SVD
【发布时间】：2020-07-12 08:27:00
【问题描述】：

我正在尝试实现上双对角矩阵 B 的分而治之 SVD，但我的代码不起作用。错误是：

"无法执行赋值，因为左边的大小是 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

【问题讨论】：

这个B1 = B(1:k-1,1:k); 应该是B1 = B(1:k,1:k); 吗？此外，您的代码不正确，例如对于B = [-1]（S 矩阵应该只包含正值，您需要使U 或V 成为B 的符号，并且S 是1x1 矩阵的B 的绝对值)
你用的B尺寸是多少？我无法复制您的错误
@Yuval Harpaz，我使用了 6×6 矩阵
现在我得到了错误，@chtz 关于 B1 = B(1:k,1:k) 的评论呢？如果您从代码中删除所有 -1，则不会出现错误。如果你不这样做，你永远不会使用 B(k,k)

标签： matlab numeric numerical-methods numerical-computing

【解决方案1】：

由于数学知识有限，您需要更多地帮助我——因为我无法以数学方式判断该方法是否正确（没有给出理论；））。无论如何，我什至无法重现错误，例如这个矩阵，MathWorks 用它来说明他们的LU matrix factorization

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

【讨论】：