I am implementing the svd decomposition using givens rotation as follows
http://reference.kfupm.edu.sa/content/b/i/the_bidiagonal_singular_value_decomposit_1343325.pdf
where:
A input matrix
U,S,V output matrices
B bidiagonal matrix
d diagonal of the bidiagonal matrix
e superdiagonal of the bidiagonal matrix
Here is the matlab code:
function [d,e, U, V]=QRZeroShift(d, e, U, V )
%Compute QR zero shift
%Algorithm by Demel and Kahan, Accurate singular values of bidiagonal %matrices, 1990, pp. 14
n = length(d);
c_old = 1;
c = 1;
for i=1:n-1
%Compute Givens rotation
[c, s, r] = givens (c *d(i), e(i));
%Update items of V matrix applying Givens rotations
[V(:, i), V(:, i + 1)] = update (c, s, V(:, i), V (:, i + 1))
if (i ~=1 )
e(i-1) = r*s_old;
end
%Compute Givens rotation
[c_old, s_old, d(i)] = givens (c_old*r,d(i+1)*s);
%Update items of U matrix applying Givens rotations
%Possible error in article: update (c, s,...) found
[U(:, i), U(:, i + 1)] = update (c_old, s_old, U(:, i), U (:, i + 1) );
end
%Actualize diagonal and superdiagonal
h = c*d(n);
e(n-1) = h*s_old;
d(n) = h*c_old;
Givens rotations:
function [c,s,r]=givens(f,g)
%Algorithm by Demel and Kahan, Accurate singular values of bidiagonal %matrices, 1990, pp. 14
if (f==0)
c=0; s=1; r=g;
elseif (abs(f) > abs (g))
t = g/f; t1 = (1 + t*t)^0.5;
c = 1/t1;s = t*c;r = f*t1;
else
t = f /g; t1 = (1 +t*t)^0.5;
s = 1/t1; c = t * s; r = g * t1;
end
Update function:
function [v1, v2 ] = update ( cs, sn, v1, v2 )
n = size(v1);
%Algorithm by Demel and Kahan
%Update procedure, replace v1->cs * t + sn * v2 (1, i)
%and v2-> -sn * t +cs * v2 (1, i)
for i=1:n
t = v1( i, 1);
v1(i, 1) = cs * t + sn * v2 (i, 1);
v2(i, 1) = -sn * t +cs * v2 (i, 1);
end
The problem is, that algorithm produces good S and V matrix, but wrong U matrix.
Example:
A = [1 2 3 4; 1 4 5 6; 30 5 7 0; 1 -2 5 6]
Biadiagonal decompoition:
U =
-0.0333 -0.4447 -0.0203 0.8949
-0.0333 -0.8089 -0.4174 -0.4127
-0.9983 0.0545 -0.0157 -0.0104
-0.0333 -0.3807 0.9084 -0.1698
d =
-30.0500
7.9476
-3.5177
-0.5285
e =
0
9.0343
8.5741
3.7482
and
V =
1.0000 0 0 0
0 -0.5673 0.1653 -0.8068
0 -0.8214 -0.0434 0.5687
0 -0.0589 -0.9853 -0.1604
After 100 iterations
for i = 1: 100
[ d, e, U, V] = QRZeroShift( d, e, U, V, 1);
end
my implementation bring the following results
d =
-31.4748
11.7767
4.3094
0.2780
U =
-0.0704 -0.4233 -0.0401 -0.9024
-0.0999 -0.6664 -0.6511 0.3493
-0.9902 0.1381 0.0169 0.0117
-0.0675 -0.5980 0.7578 0.2522
V =
0.9514 0.2085 0.1332 0.1836
0.1702 -0.1378 -0.9756 -0.0179
0.2535 -0.5628 0.1379 -0.7746
0.0409 -0.7879 0.1073 0.6050
but the correct results are:
U =
-0.0704 -0.4167 -0.2262 -0.8776
-0.0999 -0.6695 -0.5652 0.4716
-0.9902 0.1380 0.0191 0.0090
-0.0675 -0.5993 0.7931 0.0855
S=
31.4748 0 0 0
0 11.7767 0 0
0 0 4.3094 0
0 0 0 0.2780
V =
-0.9514 0.2085 0.1332 -0.1836
-0.1702 -0.1378 -0.9756 0.0179
-0.2535 -0.5628 0.1379 0.7746
-0.0409 -0.7879 0.1073 -0.6050
Matrices S, V are correct, but U is wrong... Could anybody help me? Thanks.