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.