# SVD, givens rotation, code check

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]

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.

-