SVD Matlab Implementation
Asked Answered
D

1

5

I tried to write matlab code that would decompose a matrix to its SVD form.

"Theory":

To get U, I found the eigenvectors of AA', and to get V, I found the eigenvectors of A'A. Finally, Sigma is a matrix of the same dimension as A, with the root of the eigenvalues on the diagonal in an ordered sequence.

However, it doesn't seem to work properly.

A=[2 4 1 3; 0 0 2 1];

% Get U, V
[aatVecs, aatVals] = eig(A*A');
[~, aatPermutation] = sort(sum(aatVals), 'descend');
U = aatVecs(:, aatPermutation);

[ataVecs, ataVals] = eig(A'*A);
[~, ataPermutation] = sort(sum(ataVals), 'descend');
V = ataVecs(:, ataPermutation);

% Get Sigma
singularValues = sum(aatVals(:, aatPermutation)).^0.5;
sigma=zeros(size(A));
for i=1:nnz(singularValues)
    sigma(i, i) = singularValues(i);
end

A
U*sigma*V'

U * sigma * V' seem to be returned with a factor of -1:

ans =

-2.0000   -4.0000   -1.0000   -3.0000
0.0000    0.0000   -2.0000   -1.0000

What's the mistake in the code or "theory" that led to it?

Deca answered 27/2, 2016 at 2:47 Comment(1)
I wrote about this last May: Calculating SVD by hand.Hobo
S
7

The eigenvectors are not unique (because Av==λv by definition, any w with μw==v and μ~=0 is also an eigenvector). It so happens that the eigenvectors returned by eig don't match up in the right way for SVD (even though they are normalized).

However, we can construct U once we have V, which we will find as eigenvectors of A'*A as in your algorithm. Once you have found V as sorted eigenvectors, you have to find U to match. Since V is orthogonal, A*V == U*sigma. So we can set

U = [A*V(:,1)./singularValues(1) A*V(:,2)./singularValues(2)];

and indeed, A == U*sigma*V', and in particular the U found here is exactly the negative of U found with your algorithm.

Sleepwalk answered 27/2, 2016 at 7:42 Comment(3)
While eigenvectors are not unique, the normalized eigenvectors should be, no?Deca
@AnaM No, because norm(v) == norm(-v), and with n eigenvectors you have 2^n combinations of orthonormal vectors.Sleepwalk
Wouldn't some convention to have first element as positive fix this issue?Barrow

© 2022 - 2024 — McMap. All rights reserved.