SVD computation using eig function

I'm trying to find the Singular Value Decomposition of a rectangular matrix by calculating the eigenvalues and eigenvectors of and $A^TA$. Here's what I have

clear all; close all; clc;
% Control rng
rng('default')
s = rng
A = randi([0, 20], [7,4]); % generate random 7x4 matrix A 
[U1,S1,V1] = svd(A); % Matlab's svd
% My svd
[U_tilde,L1] = eig(A*A');
[V_tilde,L2] = eig(A'*A);
% Sort eigenvalues and the corresponding eigenvectors
[D1,ind1] = sort(diag(L1),'descend');
U_tilde = U_tilde(:,ind1);
[D2,ind2] = sort(diag(L2),'descend');
V_tilde = V_tilde(:,ind2);
D = L2(ind2,ind2);
[m,n]=size(A);
L = zeros(m,n);
L(1:n,1:n) = D;

But when I compute (U_tilde)*(L.^0.5)*(V_tilde') I don't get the original matrix A.

>> U_tilde*L.^0.5*V_tilde'

ans =

    0.7203   13.7281    5.4162   21.1586
    4.8700   13.3098   26.8486   11.5016
   19.3466    9.7884   17.7988    6.4885
    2.5526   27.0760    7.0479   11.7346
   13.7950   17.3916   16.3837   16.9942
   22.6812   14.4418    9.8290   14.4702
   10.1597   11.1548    5.3215   10.1012
   
>> A

A =

    17    11    16     0
    19    20     2    17
     2    20     8    19
    19     3    19    14
    13    20    16    15
     2    20    20    15
     5    10    13     8

V1 and V_tilde are equal up to the sign, U1 and U_tilde are equal up to the 4-th column. S1 and L.^0.5 are equal. What's wrong with my code?

Answer 1

[Christine Tobler]

Try the formula the other way around, using U_tilde, V_tilde and A to compute D:

>> U_tilde'*A*V_tilde
ans =
   70.0118   -0.0000   -0.0000   -0.0000
    0.0000  -23.2868    0.0000   -0.0000
   -0.0000   -0.0000  -18.3412    0.0000
   -0.0000    0.0000   -0.0000  -11.5184
    0.0000   -0.0000   -0.0000    0.0000
    0.0000   -0.0000    0.0000    0.0000
    0.0000   -0.0000    0.0000   -0.0000

>> sqrt(diag(D))
ans =
   70.0118
   23.2868
   18.3412
   11.5184

You need to adjust the signs of the columns of either U or V.

Note that there are other problems with this approach: If A has two identical singular values, there will be two vectors U that span the left singular subspace of that singular value, and two vectors V which span the right singular subspace. But these will likely not match each other.

[Duc Anh Le]

How do we know which columns of U or V we should change signs in the general case?

[Christine Tobler]

I would really not recommend this algorithm for computing the SVD, computing the SVD directly is more accurate and also faster. Are you planning to use this in practice, or is it just a theoretical consideration?

You could compute diag(U_tilde'*A*V_tilde), find all negative values, and then flip the signs of the corresponding columns of either U_tilde or V_tilde, it doesn't matter which.

[Duc Anh Le]

Thanks for the answer, I'm studying the Divide&Conquer Algorithm for SVD and I need to implement it in MATLAB, but I can't seem to find any resources on that.

[Christine Tobler]

There's a paper about the algorithm: "A divide-and-conquer algorithm for the bidiagonal SVD" M. Gu, S.C.Eisenstat, SIAM Journal on Matrix Analysis and Applications, 1995. If you search for this on scholar.google.com, there's a link to the PDF there.