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

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.