我试图弄清楚如何最好地执行稀疏矩阵和稀疏向量的元素加法(和减法)。我在 SO 上找到了 this trick:

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

谢谢

最佳答案

使用 10x10 稀疏垫和 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 格式的变体,在更改稀疏结构时应该更有效:
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

总和有 38 个非零项,高于原始 mat 中的 10 个。它添加了 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 比 2x 好
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-7 倍。

对于非常大的 mat ,内存问题可能会影响密集性能,但迭代稀疏解决方案也不会发光。

我可以通过更有效地索引 addvec2 来从 Ml 中榨取更多性能。 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

关于python - 具有广播的稀疏 scipy 矩阵向量的元素相加,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/30741461/

10-13 02:44