# SVD computation using eig function

19 views (last 30 days)
Duc Anh Le on 20 Apr 2020
Commented: Christine Tobler on 22 Apr 2020
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?

Christine Tobler on 20 Apr 2020
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.

Show 1 older comment
Christine Tobler on 21 Apr 2020
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 on 21 Apr 2020
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 on 22 Apr 2020
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.