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

### Rohit Maheshwari (view profile)

on 14 Sep 2018
Latest activity Commented on by Rohit Maheshwari

on 16 Sep 2018

### John D'Errico (view profile)

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);

### Release

R2015a ### John D'Errico (view profile)

on 14 Sep 2018
Edited by John D'Errico

### John D'Errico (view profile)

on 14 Sep 2018

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.

Show 1 older comment
Rohit Maheshwari

### Rohit Maheshwari (view profile)

on 16 Sep 2018
Could you tell me how you came up with the logic to limit the number of iterations? Ideally, we should've stopped when matrix A has been diagonalized but on MATLAB this leads to an infinite loop.
John D'Errico

### John D'Errico (view profile)

on 16 Sep 2018
OK. Lets first agree on one thing. If A is already a diagonal matrix, then nothing is changed by the pair of QR factorizations. Agreed?
Next, you should recognize that a characteristic of a QR (due to the Householder transformations a QR is built from), is it tends to move mass onto the diagonal.
A = rand(3,2)
A =
0.86031 0.42864
0.5009 0.26735
0.99476 0.47295
>> [Q,R] = qr(A)
Q =
-0.61131 -0.21856 -0.76062
-0.35592 -0.7825 0.5109
-0.70684 0.58304 0.40055
R =
-1.4073 -0.69149
0 -0.027132
0 0
>> [Q,R] = qr(R')
Q =
-0.89751 -0.44099
-0.44099 0.89751
R =
1.568 0.011965 0
0 -0.024351 0
>> [Q,R] = qr(R')
Q =
-0.99997 -0.0076303 0
-0.0076303 0.99997 0
0 0 1
R =
-1.5681 0.00018581
0 -0.024351
0 0
>> [Q,R] = qr(R')
Q =
-1 0.0001185
0.0001185 1
R =
1.5681 -2.8855e-06 0
0 -0.024351 0
>> [Q,R] = qr(R')
Q =
-1 1.8401e-06 0
1.8401e-06 1 0
0 0 1
R =
-1.5681 -4.4808e-08
0 -0.024351
0 0
So what you see, is the stuff that goes into the upper triangle at each iteration gets smaller. It NEVER goes away. So you need a stopping point. A good point to stop is when that stuff is small, and sufficiently small that it is essentially insignificant. At the same time, you should see that each Q matrix is approaching an identity matrix.
Does this simple SVD have a problem? Well, yes. It can be made to fail.
A = [diag([3 2 2 1]);zeros(2,4)];
A = orth(rand(6))*A*orth(rand(4));
svd(A)
ans =
3
2
2
1
So svd has no problems here. But if you try the iterative QR, it will grind to a halt, with the off-diagonal stuff not converging to zero. Hey, I did say it had some issues.
As I recall, the simple algorithm will converge as long as the singular values are distinct. If a pair of singular values is close to each other, then convergence will be slower than otherwise. The trick would then be how to accelerate convergence for those problem matrices. (Something that arguably goes way beyond the scope of Answers comments, and into a class on numerical linear algebra.) The simple solution is to just use the carefully written SVD that we already have.
Rohit Maheshwari

### Rohit Maheshwari (view profile)

on 16 Sep 2018
Once again, thank you for the clarification.