用廣播元素添加稀疏的scipy矩陣矢量

Question

我想弄清楚如何最好地執行稀疏矩陣和稀疏矢量的元素加法（和減法）。我發現this trick上SO：

mat = sp.csc_matrix([[1,0,0],[0,1,0],[0,0,1]]) 
vec = sp.csr_matrix([[1,2,1]]) 
mat.data += np.repeat(vec.toarray()[0], np.diff(mat.indptr)) 

但不幸的是它僅更新非零值：

print(mat.todense()) 
[[2 0 0] 
[0 3 0] 
[0 0 2]] 

的SO線程上的實際接受的答案：

def sum(X,v): 
    rows, cols = X.shape 
    row_start_stop = as_strided(X.indptr, shape=(rows, 2), 
          strides=2*X.indptr.strides) 
    for row, (start, stop) in enumerate(row_start_stop): 
     data = X.data[start:stop] 
     data -= v[row] 

sum(mat,vec.A[0]) 

做同樣的事情。我現在不幸沒有想法，所以我希望你能幫我找出解決這個問題的最好方法。

編輯： 我期待它做一樣的這種密集的版本會做：

np.eye(3) + np.asarray([[1,2,1]]) 
array([[ 2., 2., 1.], 
     [ 1., 3., 1.], 
     [ 1., 2., 2.]]) 

感謝

Answer 1

所以你的原始矩陣是稀疏的，載體是稀疏，但在由此產生的矩陣與你的向量中的非零座標對應的列將是密集的。

所以我們不妨兌現這些列的稠密矩陣

def addvec(mat,vec): 
    for i in vec.nonzero()[1]: 
     mat[:,i] = sp.csc_matrix(mat[:,i].todense() + vec[0,i]) 
    return mat 

Answer 2

與10×10稀疏墊和VEC一些測試：

In [375]: mat=sparse.rand(10,10,.1) 
In [376]: mat 
Out[376]: 
<10x10 sparse matrix of type '<class 'numpy.float64'>' 
    with 10 stored elements in COOrdinate format> 

In [377]: mat.A 
Out[377]: 
array([[ 0.  , 0.  , 0.  , 0.  , 0.  , 
     0.  , 0.  , 0.  , 0.  , 0.  ], 
     [ 0.  , 0.  , 0.  , 0.  , 0.  , 
     0.15568621, 0.59916335, 0.  , 0.  , 0.  ], 
     ... 
     [ 0.  , 0.  , 0.15552687, 0.  , 0.  , 
     0.47483064, 0.  , 0.  , 0.  , 0.  ]]) 

In [378]: vec=sparse.coo_matrix([0,1,0,2,0,0,0,3,0,0]).tocsr() 
<1x10 sparse matrix of type '<class 'numpy.int32'>' 
    with 3 stored elements in Compressed Sparse Row format> 

maxymoo的解決方案：

def addvec(mat,vec): 
    Mc = mat.tocsc() 
    for i in vec.nonzero()[1]: 
     Mc[:,i]=sparse.csc_matrix(Mc[:,i].todense()+vec[0,i]) 
    return Mc  

和變異使用lil格式，當c時應該更有效掛稀疏結構：

def addvec2(mat,vec): 
    Ml=mat.tolil() 
    vec=vec.tocoo()            
    for i,v in zip(vec.col, vec.data): 
     Ml[:,i]=sparse.coo_matrix(Ml[:,i].A+v) 
    return Ml 

的Sumation公司有38個非零項，從10中的原始mat。它添加了來自vec的3列。這是稀疏性的一個重大變化。

In [382]: addvec(mat,vec) 
Out[382]: 
<10x10 sparse matrix of type '<class 'numpy.float64'>' 
    with 38 stored elements in Compressed Sparse Column format> 

In [383]: _.A 
Out[383]: 
array([[ 0.  , 1.  , 0.  , 2.  , 0.  , 
     0.  , 0.  , 3.  , 0.  , 0.  ], 
     [ 0.  , 1.  , 0.  , 2.  , 0.  , 
     0.15568621, 0.59916335, 3.  , 0.  , 0.  ], 
     ... 
     [ 0.  , 1.  , 0.15552687, 2.  , 0.  , 
     0.47483064, 0.  , 3.  , 0.  , 0.  ]]) 

與addvec2相同輸出：

In [384]: addvec2(mat,vec) 
Out[384]: 
<10x10 sparse matrix of type '<class 'numpy.float64'>' 
    with 38 stored elements in LInked List format> 

而在定時，addvec2確實優於2×

In [385]: timeit addvec(mat,vec) 
100 loops, best of 3: 6.51 ms per loop 

In [386]: timeit addvec2(mat,vec) 
100 loops, best of 3: 2.54 ms per loop 

和緻密當量：

In [388]: sparse.coo_matrix(mat+vec.A) 
Out[388]: 
<10x10 sparse matrix of type '<class 'numpy.float64'>' 
    with 38 stored elements in COOrdinate format> 

In [389]: timeit sparse.coo_matrix(mat+vec.A) 
1000 loops, best of 3: 716 µs per loop 

In [390]: timeit sparse.coo_matrix(mat.A+vec.A) 
1000 loops, best of 3: 338 µs per loop 

一個版本威力節省臨時密集矩陣空間，運行在同一時間：

In [393]: timeit temp=mat.A; temp+=vec.A; sparse.coo_matrix(temp) 
1000 loops, best of 3: 334 µs per loop 

那麼密集的版本確實5-7x比我稀疏的版本更好。

對於一個非常大的mat，內存問題可能會咀嚼密集的性能，但迭代稀疏解決方案也不會發光。

我可以通過更有效地索引Ml來從addvec2中獲得更多性能。 Ml.data[3],Ml.rows[3]比Ml[3,:]或Ml[:,3]快得多。

def addvec3(mat,vec): 
    Mtl=mat.T.tolil() 
    vec=vec.tocoo() 
    n = mat.shape[0] 
    for i,v in zip(vec.col, vec.data): 
     t = np.zeros((n,))+v 
     t[Mtl.rows[i]] += Mtl.data[i] 
     t = sparse.coo_matrix(t) 
     Mtl.rows[i] = t.col 
     Mtl.data[i] = t.data 
    return Mtl.T 

In [468]: timeit addvec3(mat,vec) 
1000 loops, best of 3: 1.8 ms per loop 

適度的改善，但沒有我希望的那麼多。並擠多一點：

def addvec3(mat,vec): 
    Mtl = mat.T.tolil() 
    vec = vec.tocoo(); 
    t0 = np.zeros((mat.shape[0],)) 
    r0 = np.arange(mat.shape[0]) 
    for i,v in zip(vec.col, vec.data): 
     t = t0+v 
     t[Mtl.rows[i]] += Mtl.data[i] 
     Mtl.rows[i] = r0 
     Mtl.data[i] = t 
    return Mtl.T 

In [531]: timeit mm=addvec3(mat,vec) 
1000 loops, best of 3: 1.37 ms per loop