Numpy svd和Scipy.sparse svds - python

我正在为Python中的稀疏未定系统实现一个求解器(讨论了here),并且尝试使用scipy在SciPy cookbook中重建使用标准numpy svd函数(numpy.linalg.svd)的nullspace函数。 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是一个csc_matrix

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]]

而以下

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的原因。

参考方案

您发布的问题中的输出对我来说看起来不错。在numpy调用中,您正在计算每个奇异值,在scipy代码中,您仅计算了前k个奇异值,它们与numpy输出中的前k个匹配。

稀疏的前k个svd不会让您计算每个奇异值,因为如果您想这样做,则可以使用完整的svd函数。

下面,我提供了一些代码供您自己检查。需要注意的是,虽然numpy和scipy完整svds都可以很好地重新创建原始矩阵,但前k个svd却不能。这是因为您要丢弃数据。通常情况下,如果您可以保持足够的距离就可以了。问题在于,如果将SVD与前k个一起使用,它可以用作原始矩阵的低秩近似,而不是替代。

为了清楚起见,我的经验来自为原作者A Sparse Plus Low-Rank Exponential Language Model for Limited Resource Scenarios实现本文的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)

用Python计算稀疏矩阵的Cholesky分解 - python

我正在尝试实现Reinsch's Algorithm(pp 4)。由于工作矩阵是稀疏的,所以我使用的是scipy.sparse模块,但是正如您所看到的,Reinsch的算法需要稀疏矩阵的Cholesky分解(我们称其为my_matrix)才能求解某些系统,但是我不能找到与此有关的任何东西。当然,在同一算法中,我可以使用scipy.sparse.li…

Python:如何停止多线程的numpy? - python

我知道这似乎是一个荒谬的问题,但是我必须在与部门中其他人共享的计算服务器上定期运行作业,当我开始10个作业时,我真的希望它只占用10个核心而不是更多;我不在乎每次运行一个内核所需的时间是否更长:我只是不想让它侵犯其他人的领土,这将需要我放弃工作等等。我只想拥有10个核心,仅此而已。更具体地说,我在基于Python 2.7.3和numpy 1.6.1的Redh…

numpy.savetxt“元组索引超出范围”? - python

我试图在文本文件中写几行,这是我使用的代码:import numpy as np # Generate some test data data = np.arange(0.0,1000.0,50.0) with file('test.txt', 'w') as outfile: outfile.write('…

numpy loadtxt单行/行作为列表 - python

我只有一个数据文件,例如: 1.2 2.1 3.2 我使用numpy版本1.3.0 loadtxt加载它 a,b,c = loadtxt("data.dat", usecols(0,1,2), unpack=True) 输出是浮点数而不是数组 a = 1.2 我希望它将是: a = array([1.2]) 如果我读取了多行文件,则该文件…

python / numpy中的多线程blas - python

我正在尝试在Python中实现大量矩阵矩阵乘法。最初,我假设NumPy将自动使用我的线程BLAS库,因为我是根据这些库构建它的。但是,当我查看top或其他内容时,似乎代码根本不使用线程。有什么想法是错误的,或者我可以做些什么来轻松使用BLAS性能? 参考方案 并非所有的Nu​​mPy都使用BLAS,只有某些功能-特别是dot(),vdot()和innerpr…