MATLAB: Where am I going wrong in this code to find SVD of a matrix

MATLABorthogonal matrixsingular value decomposition

I'm working on my own code to find the Singular Value Decomposition of a square matrix. The answer I'm getting doesn't match the actual answer. What changes do I need to make to my code?
A=input('Enter matrix A:\n');
C=A*A';D=A'*A;
[v1,e1]=eig(C);
[v2,e2]=eig(D);
n=max(size(C));
j=n;
for i=1:n/2;
    t=v1(:,i);
    v1(:,i)=v1(:,j);
    v1(:,j)=t;
    j=j-1;
end
U=orth(v1);
display('Matrix U is:');
display(U);
j=n;
for i=1:n/2;
    t=v2(:,i);
    v2(:,i)=v2(:,j);
    v2(:,j)=t;
    j=j-1;
end
V=orth(v2);
display('Matrix V is:');
display(V);

Best Answer

Why do I find this task a little silly? At least, the way you are accomplishing it is driving the wrong way down a one-way street.
You call the MATLAB function orth TWICE in the code, in an attempt to compute the SVD of a matrix. Do you realize that orth itself calls svd?
So, effectively, you are using svd to compute the svd. But worse, you are using svd TWICE to compute the svd. Not only that, but you also call eig twice. Are you trying to save time? ;-)
Sorry, but I gotta laugh. No, there is nothing overtly silly about wanting to understand the SVD, even in writing code to compute it from scratch in MATLAB. But using svd to compute the svd?
Were I to write an SVD code from scratch, I might use a basic scheme, perhaps Householder transformations to bidiagonalize the matrix. Then further rotations to kill off the off-diagonal elements. I recall that from many years ago, when I did write my own svd. I'd imagine that some light reading might show there are now better ways to accomplish the task.
Now, what did you do wrong? You say that you don't get the correct result. Why do you think it is wrong? Is it possible that the singular vectors you found are the same, yet they have sign changes? Or, perhaps, is your matrix numerically singular? In that case, the vectors spanning the null space are not unique. Worse, when you compute C and D
C=A*A';D=A'*A;
you kill off the small singular values of A, pushing them into the null space. So any singular values less than sqrt(eps) times the maximum singular value will now no longer be computable.
So there are MANY problems in what you are doing. You CAN do it right, but that means you need to start from basics, and not try to use high level codes to recompute the svd. One idea might be to use QR, iteratively. The following code is just a quick hack, but it works reasonably on simple matrices.
function [U,S,V] = mysvd(A)
[n,m] = size(A);
U = eye(n);
V = eye(m);
tol = max(abs(A(:)))*1.e-15;
Arem = inf;
while Arem > tol
  [Qu,Ru] = qr(A);
  U = U*Qu;
    [Qv,Rv] = qr(Ru');
    V = V*Qv;
    A = Rv';
    Arem = norm(tril(A,-1),inf);
  end
U = U(:,1:m);
S = triu(A(1:m,:));
% correct any negative singular values
V = V.*sign(diag(S)).';
S = abs(S);
Does it work? A test is always good.
A = rand(5,3);
[U,S,V] = mysvd(A)
U =
     -0.25721      0.41512    0.0046487
     -0.35691    -0.002447      0.91251
     -0.71507    -0.071517     -0.40011
     -0.37533     -0.72126    -0.023075
     -0.39275      0.54986    -0.081758
S =
       1.8273            0            0
            0       0.7709            0
            0            0      0.57832
V =
     -0.65496     -0.72374     -0.21734
     -0.67233       0.6894     -0.26963
     -0.34497     0.030467      0.93812
>> U*S*V' - A
ans =
   3.1919e-16  -1.1102e-15  -1.9429e-16
   6.1062e-16  -6.6613e-16  -2.2204e-16
  -8.8818e-16  -8.8818e-16   1.9429e-16
  -7.7716e-16  -1.3878e-17  -5.5511e-17
   4.1633e-16  -1.4433e-15   8.3267e-17
Compare the results to svd.
[U0,S0,V0] = svd(A,0)
U0 =
     -0.25721     -0.41512    0.0046487
     -0.35691     0.002447      0.91251
     -0.71507     0.071517     -0.40011
     -0.37533      0.72126    -0.023075
     -0.39275     -0.54986    -0.081758
S0 =
       1.8273            0            0
            0       0.7709            0
            0            0      0.57832
V0 =
     -0.65496      0.72374     -0.21734
     -0.67233      -0.6894     -0.26963
     -0.34497    -0.030467      0.93812
Essentially correct to within sign differences in the singular vectors, with no need to resort to anything more sophisticated than repeated calls to QR. And if you wanted to do so, you could easily replace the QR calls with Householder transformations. Even better, the code I wrote should be numerically stable, since there was no need to create A'*A, etc.
The code I wrote fails to run if A has more columns than rows, but even that is trivially repaired. I have a funny feeling it might also be screwed up with matrices where the singular values are not distinct. But hey, that too can be repaired with some effort.
Related Question