Divide and Conquer SVD in MATLAB

Question

I'm trying to implement Divide and Conquer SVD of an upper bidiagonal matrix B, but my code is not working. The error is:

"Unable to perform assignment because the size of the left side is 3-by-3 and the size of the right side is 2-by-2.
V_bar(1:k,1:k) = V1;"

Can somebody help me fix it? Thanks.

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

Should this `B1 = B(1:k-1,1:k);` be `B1 = B(1:k,1:k);`? Also, your code is incorrect, e.g. for `B = [-1]` (the `S` matrix should contain only positive values, you need to make either `U` or `V` be the sign of `B`, and `S` the absolute value of `B` for 1x1 matrices) — chtz, Mar 31 '20 at 18:26
what is the size of B you use? I was unable to replicate your error — Yuval Harpaz, Apr 05 '20 at 16:14
now I get the error, what about @chtz's comment about B1 = B(1:k,1:k)? if you remove all the -1 from your code you get no error. ANd if you don't, you never use B(k,k) — Yuval Harpaz, Apr 06 '20 at 05:40

score 0 · Answer 1 · answered Apr 05 '20 at 07:49

With limited mathematical knowledge, you need to help me a bit more -- as I cannot judge if the approach is correct in a mathematical way (with no theory given;) ). Anyway, I couldn't even reproduce the error e.g with this matrix, which The MathWorks use to illustrate their LU matrix factorization

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

So I tried to structure your code a bit and gave some hints. Extend this to make your code clearer for those people (like me) who are not too familiar with matrix decomposition.

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

Divide and Conquer SVD in MATLAB

1 Answers1