为什么Octave、R、Numpy和LAPACK在同一个矩阵上产生不同的SVD结果?

9
我将使用Octave和R来计算一个简单矩阵的SVD,但是得到了两个不同的答案!以下是代码:

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

我对Octave和R都很陌生,所以我的第一个问题是:我是否做得对?如果是这样,哪一个是“正确”的?它们两个都正确吗?我也尝试了在numpy中调用LAPACK函数dgesdd和dgesvd。Numpy给我一个类似于Octave的答案,而调用LAPACK函数则给我一个类似于R的答案。

谢谢!

2个回答

8
在SVD分解$A=UDV^T$中,只有$D$是唯一的(可以重新排序)。很容易看出$cU$和$\frac{1}{c}V$会给出相同的分解。因此,不同的算法可能会得出不同的结果,重要的是所有算法的$D$必须相同。

5
很抱歉,这里没有TeX支持。 - chl
3
只有当你采用对角线元素按降序排列的约定时,矩阵 D 才能唯一确定。 - Sven Marnach
1
另一个问题是数值精度,即在1e-16数量级的数字实际上等同于零。 - IRTFM
@chl,真遗憾。有什么想法吗? - mpiktas
2
@mpiktas 不是的。啊,是一个故意的特性,鼓励你在{maths|stats}.SE上做出贡献 :-) - chl
显示剩余2条评论

2
实际上,对于唯一的奇异值,U和V也是唯一的。你的问题在于奇异值2和3重复了。因此,奇异值1(即3.4)是唯一的,因此U和V的第一列在两个答案中是相同的。
此外,即使第二列和第三列不是唯一的,它们应该在两个答案中处于相同的线性子空间中。这意味着如果U1由第一个U的第二列和第三列组成,U2由第二个U的第二列和第三列组成,则U1'* U2应该是一个完整的秩为2x2矩阵,其列都具有单位(欧几里得)幅度。

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接