Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

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.

share|improve this question
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.