Numpy svd vs Scipy.sparse svds

Question

I was working on implementing a solver for sparse undetermined systems in Python (discussed here) and I was trying to rebuild the nullspace function that uses the standard numpy svd function (numpy.linalg.svd) in the SciPy cookbook using the scipy.sparse version of svd (scipy.sparse.linalg.svds) but it outputs different left and right singular vectors for the examples I ran - including the matrices:

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

Why do these two solvers produce two different svd outputs for the matrices above? What can I do to ensure the same output?

Edit

Here's an example: table is a csc_matrix such that

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

So, the following code outputs

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

Note that there are quite a few values that overlap between the two solutions. Also, the scipy.sparse.linalg.svds function doesn't allow k to be greater than or equal to min(table.shape), which is why I chose k=2.

Answer 1

The output in the question you posted looks fine to me. In the numpy call you are calculating every singular value and in the scipy code you are calculating just the top k singular values, and they match the top k from the numpy output.

The sparse top k svd won't let you calculate every singular value because if you wanted to do that, then you could just use the full svd function.

Below I have included code for you to check this out yourself. The caveat is that while the numpy and scipy full svds can both recreate the original matrix well enough, the top k svd cannot. This is because you are throwing away data. Normally this is fine given that you are okay with being close enough. The issue is that SVD if used with top k can be used as a low rank approximation of the original matrix, not as a replacement.

For clarity, my experience on this comes from implementing a python, parallel version of this paper for the original author, A Sparse Plus Low-Rank Exponential Language Model for Limited Resource Scenarios.

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)