Numpy svd和Scipy.sparse svds的比较

Question

我正在实现Python中稀疏超定系统的求解器（在此讨论：这里），并尝试重建使用标准numpy svd函数（numpy.linalg.svd）的nullspace函数在SciPycookbook中的实现，但是我发现它使用scipy.sparse版本的svd（scipy.sparse.linalg.svds）为我运行的示例输出不同的左奇异向量和右奇异向量，包括以下矩阵：

[[1,1,0,0,0],[0,0,1,1,0],[1,1,1,1,1]]  
[[1,0,1],[1,1,0],[0,1,1]]

这两个求解器为何会针对上述矩阵产生不同的SVD输出？我该如何确保相同的输出结果？

编辑

以下是一个示例： table 是一个，具体如下：

table.todense()  = matrix([[1,1,0,0,0],[0,0,1,1,0],[1,1,1,1,1]],dtype=int64)

因此，以下代码会输出：

numpy.linalg.svd(table.todense()) =  
[[ -3.64512933e-01   7.07106781e-01  -6.05912800e-01]  
[ -3.64512933e-01  -7.07106781e-01  -6.05912800e-01]  
[ -8.56890100e-01   2.32635116e-16   5.15499134e-01]]  
-----------------------------------------------------
[ 2.58873755  1.41421356  0.54629468]
-----------------------------------------------------
[[ -4.7181e-01 -4.7181e-01 -4.7181e-01 -4.7181e-01 -3.3101e-01]
[5e-01   5e-01  -5e-01  -5e-01 6.16450329e-17]
[-1.655e-01  -1.655e-01  -1.655e-01  -1.655e-01  9.436e-01]
[5e-01  -5e-01  -5e-01   5e-01 -1.77302319e-16]
[-5e-01  5e-01  -5e-01   5e-01 2.22044605e-16]]

And the following

scipy.sparse.linalg.svds(table,k=2)=  
[[  7.07106781e-01,  -3.64512933e-01],
[ -7.07106781e-01,  -3.64512933e-01],
[  2.73756255e-18,  -8.56890100e-01]]
-------------------------------------
[ 1.41421356,  2.58873755]
-------------------------------------
[[  5e-01,   5e-01,  -5e-01, -5e-01,   1.93574904e-18],
[ -4.71814e-01,  -4.71814e-01,  -4.71814e-01, -4.71814e-01,  -3.31006e-01]]

请注意，两个解决方案之间存在许多重叠的值。此外，scipy.sparse.linalg.svds函数不允许k大于或等于min(table.shape)，这就是我选择k = 2的原因。

- Rushabh Mehta

---

请问您能否提供一组输入数据、演示如何运行这些函数以及输出结果的示例呢？这将极大地改善这个问题。 - AGN Gazer

我建议你查看相关联的问题，我在那里解释得更清楚。 - Rushabh Mehta

对于这两个函数，奇异值(s)是它们承诺的主要内容。前k个最大的s值匹配。SVD并不唯一。https://en.wikipedia.org/wiki/Singular-value_decomposition#Calculating_the_SVD - hpaulj

计算截断SVD存在哪些快速算法看起来是一个有用的答案。它讨论了截断SVD，并提到了ARPACK实现，svds也引用了它。MATLAB具有svd和svds两种方法。此外，http://fa.bianp.net/blog/2012/singular-value-decomposition-in-scipy/也是一个不错的参考。 - hpaulj

哦，有趣。我没想到奇异值分解不是唯一的。这肯定会对我的计划造成影响。 - Rushabh Mehta

Answer 1

你发布的问题中的输出看起来很好。在numpy调用中，你正在计算每个奇异值，在scipy代码中，你只计算前k个奇异值，并且它们与numpy输出的前k个匹配。

稀疏的前k个奇异值分解不能让你计算每个奇异值，因为如果你想这样做，那么你可以使用完整的svd函数。

下面我包括了代码，让你自己检查。但需要注意的是，虽然numpy和scipy的完整svd都可以很好地重建原始矩阵，但前k个奇异值分解却不能。这是因为你正在丢弃数据。通常这没关系，因为你可以接受足够接近。问题在于，如果使用前k个奇异值分解，则可以将其用作原始矩阵的低秩近似，而不是替代品。

为了明确起见，我的经验来自为原作者实现此论文的python并行版本，用于有限资源场景的稀疏加低秩指数语言模型。

import numpy as np                                                                                   
from scipy import linalg                                                                            
from scipy.sparse import linalg as slinalg                                                           

x = np.array([[1,1,0,0,0],[0,0,1,1,0],[1,1,1,1,1]],dtype=np.float64)                                 

npsvd = np.linalg.svd(x)                                                                             
spsvd = linalg.svd(x)                                                                                
sptop = slinalg.svds(x,k=2)                                                                          

print "np"                                                                                           
print "u: ", npsvd[0]                                                                                
print "s: ", npsvd[1]                                                                                
print "v: ", npsvd[2]                                                                                

print "\n=========================\n"                                                                

print "sp"                                                                                           
print "u: ", spsvd[0]                                                                                
print "s: ", spsvd[1]                                                                                
print "v: ", spsvd[2]                                                                                

print "\n=========================\n"                                                                

print "sp top k"                                                                                     
print "u: ", sptop[0]                                                                                
print "s: ", sptop[1]                                                                                
print "v: ", sptop[2]                                                                                

nptmp = np.zeros((npsvd[0].shape[1],npsvd[2].shape[0]))                                              
nptmp[np.diag_indices(np.min(nptmp.shape))] = npsvd[1]                                               
npreconstruct = np.dot(npsvd[0], np.dot(nptmp,npsvd[2]))                                             

print npreconstruct                                                                                  
print "np close? : ", np.allclose(npreconstruct, x)                                                  

sptmp = np.zeros((spsvd[0].shape[1],spsvd[2].shape[0]))                                              
sptmp[np.diag_indices(np.min(sptmp.shape))] = spsvd[1]                                               
spreconstruct = np.dot(spsvd[0], np.dot(sptmp,spsvd[2]))                                             

print spreconstruct                                                                                  
print "sp close? : ", np.allclose(spreconstruct, x)                                                  

sptoptmp = np.zeros((sptop[0].shape[1],sptop[2].shape[0]))                                           
sptoptmp[np.diag_indices(np.min(sptoptmp.shape))] = sptop[1]                                         
sptopreconstruct = np.dot(sptop[0], np.dot(sptoptmp,sptop[2]))                                       

print sptopreconstruct                                                                               
print "sp top close? : ", np.allclose(sptopreconstruct, x)

- David

---

这很有趣。能否通过这个截断的奇异值分解计算矩阵的零空间呢？如果不能，那么在Python中是否有任何用于稀疏矩阵的完整奇异值分解求解器？ - Rushabh Mehta

从我的经验来看，我认为你不能使用截断的奇异值分解来计算零空间，因为你需要所有的奇异值。稀疏版本和完整版本之间有两个主要区别。完整版本比稀疏版本快一个数量级（O(n^3) v O(n^4)），但在内存需求方面随着矩阵大小而扩展。稀疏版本在非零数的数量上扩展内存。在你的情况下，只要矩阵不太大，我会使用完整版本。如果需要稀疏矩阵，你也可以查看 sparsesvd 等工具。 - David

你尝试过将稀疏矩阵传递给scipy.linalg.svd吗？ - David

是的，我有...我得到了这个：此函数不支持稀疏矩阵。也许可以尝试使用scipy.sparse.linalg函数。 - Rushabh Mehta

如果它有完整的SVD，我一定愿意尝试！ - Rushabh Mehta

Answer 2

u,sigma,v = numpy.linalg.svd(A)

等于

u,sigma,v = scipy.sparse.linalg.svds(A,k=min(A.shape[0],A.shape[1]))
u = -u[:,::-1]
v = -v[::-1,:]
sigma = sigma[::-1]

- Roger Hung

---

这个回答如何解决问题？ - Rushabh Mehta

这并不是说上面的回答没有尝试回答“为什么这两个求解器会产生两个不同的SVD输出”，或者没有正确地回答这个问题 - 但是，请解释一下。请注意如何撰写一个好的答案？。 - greybeard