Any reason why Octave, R, Numpy and LAPACK yield different SVD results on the same matrix?

Question

I'm using Octave and R to compute SVD using a simple matrix and getting two different answers! The code is listed below:

R

> a<-matrix(c(1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,1,1,1), 9, 4)
> a
 [,1] [,2] [,3] [,4]
 [1,]    1    1    0    0
 [2,]    1    1    0    0
 [3,]    1    1    0    0
 [4,]    1    0    1    0
 [5,]    1    0    1    0
 [6,]    1    0    1    0
 [7,]    1    0    0    1
 [8,]    1    0    0    1
 [9,]    1    0    0    1
> a.svd <- svd(a)
> a.svd$d
[1] 3.464102e+00 1.732051e+00 1.732051e+00 1.922963e-16
> a.svd$u
       [,1]       [,2]          [,3]          [,4]
 [1,] -0.3333333  0.4714045 -1.741269e-16  7.760882e-01
 [2,] -0.3333333  0.4714045 -3.692621e-16 -1.683504e-01
 [3,] -0.3333333  0.4714045 -5.301858e-17 -6.077378e-01
 [4,] -0.3333333 -0.2357023 -4.082483e-01  6.774193e-17
 [5,] -0.3333333 -0.2357023 -4.082483e-01  6.774193e-17
 [6,] -0.3333333 -0.2357023 -4.082483e-01  6.774193e-17
 [7,] -0.3333333 -0.2357023  4.082483e-01  5.194768e-17
 [8,] -0.3333333 -0.2357023  4.082483e-01  5.194768e-17
 [9,] -0.3333333 -0.2357023  4.082483e-01  5.194768e-17
> a.svd$v
      [,1]       [,2]          [,3] [,4]
[1,] -0.8660254  0.0000000 -4.378026e-17  0.5
[2,] -0.2886751  0.8164966 -2.509507e-16 -0.5
[3,] -0.2886751 -0.4082483 -7.071068e-01 -0.5
[4,] -0.2886751 -0.4082483  7.071068e-01 -0.5

Octave

octave:32> a = [ 1, 1, 1, 1, 1, 1, 1, 1, 1; 1, 1, 1, 0, 0, 0, 0, 0, 0; 0, 0, 0, 1, 1, 1, 0, 0, 0; 0, 0, 0, 0, 0, 0, 1, 1, 1 ]
a =

  1   1   1   1   1   1   1   1   1
  1   1   1   0   0   0   0   0   0
  0   0   0   1   1   1   0   0   0
  0   0   0   0   0   0   1   1   1

octave:33> [u, s, v] = svd (a)
u =

 -8.6603e-01  -1.0628e-16   2.0817e-17  -5.0000e-01
 -2.8868e-01   5.7735e-01  -5.7735e-01   5.0000e-01
 -2.8868e-01  -7.8868e-01  -2.1132e-01   5.0000e-01
 -2.8868e-01   2.1132e-01   7.8868e-01   5.0000e-01

s =

Diagonal Matrix

 3.4641e+00            0            0            0            0
          0   1.7321e+00            0            0            0
          0            0   1.7321e+00            0            0
          0            0            0   3.7057e-16            0
          0            0            0            0            0

v =

 -3.3333e-01   3.3333e-01  -3.3333e-01   7.1375e-01
 -3.3333e-01   3.3333e-01  -3.3333e-01  -1.3472e-02
 -3.3333e-01   3.3333e-01  -3.3333e-01  -7.0027e-01
 -3.3333e-01  -4.5534e-01  -1.2201e-01  -9.0583e-17
 -3.3333e-01  -4.5534e-01  -1.2201e-01   2.0440e-17
 -3.3333e-01  -4.5534e-01  -1.2201e-01   2.0440e-17
 -3.3333e-01   1.2201e-01   4.5534e-01   8.3848e-17
 -3.3333e-01   1.2201e-01   4.5534e-01   8.3848e-17
 -3.3333e-01   1.2201e-01   4.5534e-01   8.3848e-17

I'm new to both Octave and R so my first question is am I doing this right? If so, which one is "correct"? Are they both right? I've also tried this out in numpy and calling the LAPACK functions dgesdd and dgesvd directly. Numpy give me an answer similar to Octave and calling the LAPACK functions gives me an answer similar to R.

Thanks!

Oswald

Answer 1

In SVD decomposition $A=UDV^T$ only $D$ is unique (up to reordering). It is more or less easy to see that $cU$ and $\frac{1}{c}V$ will give the same decomposition. So it is not surprising that different algorithms can give different results. What matters is that $D$ must be the same for all algorithms.

Riocard

@Riocard (+1) He, he: no TeX support here :-( – Myiasis

D is only uniquely determined if you adopt the convention that its diagonal entries are sorted in descending order. – Leasia

@Sven, yes naturaly, but if you compare two sets of numbers different ordering seems lesser issue than different values. I will edit the answer to reflect your comment. – Riocard

And the other issue is numerical accuracy, i.e., that numbers on the order of 1e-16 are effectively zero. – Gravitate

@chl, it is a pity. Any ideas why? – Riocard

@Riocard Nope. Ah yes, it is a deliberate feature to encourage your contributions on {maths|stats}.SE :-) – Myiasis

Thanks everyone for the answers! I suspected that as long A = U * D * V' all results were valid. What threw me was that since all the programs I mentioned are backed by LAPACK, I expected that I would get the same output for the same input. – Oswald

Answer 2

Actually, the U and V are also unique for unique singular values. The reason yours are not is that singular values 2 and 3 are repeated. Singular value 1 (the 3.4) is unique - and therefore columns 1 of U and V are the same in both answers.

Also, even though columns 2 and 3 are not unique, they should lie in the same linear subspace for both answers. This means that if U1 consists of columns 2 and 3 of the first U, and U2 consists of columns 2 and 3 of the second U, then U1'*U2 should be a full rank 2x2 matrix whose columns are both of unit (euclidean) magnitude.

Furtek