我正在尝试计算一个简单的点积，但保留原始矩阵中的非零值不变。玩具示例:

import numpy as np

A = np.array([[2, 1, 1, 2],
              [0, 2, 1, 0],
              [1, 0, 1, 1],
              [2, 2, 1, 0]])
B = np.array([[ 0.54331039,  0.41018682,  0.1582158 ,  0.3486124 ],
              [ 0.68804647,  0.29520239,  0.40654206,  0.20473451],
              [ 0.69857579,  0.38958572,  0.30361365,  0.32256483],
              [ 0.46195299,  0.79863505,  0.22431876,  0.59054473]])

期望的结果:

C = np.array([[ 2.        ,  1.        ,  1.        ,  2.        ],
              [ 2.07466874,  2.        ,  1.        ,  0.73203386],
              [ 1.        ,  1.5984076 ,  1.        ,  1.        ],
              [ 2.        ,  2.        ,  1.        ,  1.42925865]])

然而，实际的矩阵是稀疏的，看起来更像这样:

A = sparse.rand(250000, 1700, density=0.001, format='csr')
B = sparse.rand(1700, 1700, density=0.02, format='csr')

一个简单的方法就是使用掩码索引设置值，就像这样:

mask = A != 0
C = A.dot(B)
C[mask] = A[mask]

但是，我的原始数组非常稀疏且非常大，因此通过索引分配更改它们非常慢。转换为 lil 矩阵会有所帮助，但同样，转换本身会花费大量时间。

我猜，另一种明显的方法就是求助于迭代并跳过屏蔽值，但我不想放弃 numpy/scipy 优化数组乘法的好处。

一些说明:我实际上对某种特殊情况感兴趣，其中 B 始终是正方形，因此 A 和C 形状相同。因此，如果有一个解决方案不适用于任意数组但适合我的情况，那很好。

更新:一些尝试:

from scipy import sparse
import numpy as np

def naive(A, B):
    mask = A != 0
    out = A.dot(B).tolil()
    out[mask] = A[mask]
    return out.tocsr()


def proposed(A, B):
    Az = A == 0
    R, C = np.where(Az)
    out = A.copy()
    out[Az] = np.einsum('ij,ji->i', A[R], B[:, C])
    return out


%timeit naive(A, B)
1 loops, best of 3: 4.04 s per loop

%timeit proposed(A, B)
/usr/local/lib/python2.7/dist-packages/scipy/sparse/compressed.py:215: SparseEfficiencyWarning: Comparing a sparse matrix with 0 using == is inefficient, try using != instead.

/usr/local/lib/python2.7/dist-packages/scipy/sparse/coo.pyc in __init__(self, arg1, shape, dtype, copy)
    173                     self.shape = M.shape
    174 
--> 175                 self.row, self.col = M.nonzero()
    176                 self.data = M[self.row, self.col]
    177                 self.has_canonical_format = True

MemoryError: 

另一项更新:

无法从 Cython 中获得更多或更少的有用信息，至少在与 Python 相去甚远的情况下是这样。想法是将点积留给 scipy，并尝试尽快设置这些原始值，如下所示:

cimport cython


@cython.cdivision(True)
@cython.boundscheck(False)
@cython.wraparound(False)
cpdef coo_replace(int [:] row1, int [:] col1, float [:] data1, int[:] row2, int[:] col2, float[:] data2):
    cdef int N = row1.shape[0]
    cdef int M = row2.shape[0]
    cdef int i, j
    cdef dict d = {}

    for i in range(M):
        d[(row2[i], col2[i])] = data2[i]

    for j in range(N):
        if (row1[j], col1[j]) in d:
            data1[j] = d[(row1[j], col1[j])]

这比我之前的第一个“天真”实现(使用 .tolil())要好一些，但是按照 hpaulj 的方法，lil 可以被抛弃。也许用 std::map 之类的东西替换 python dict 会有所帮助。

最佳答案

您的 naive 代码的一个可能更干净、更快的版本:

In [57]: r,c=A.nonzero()    # this uses A.tocoo()

In [58]: C=A*B
In [59]: Cl=C.tolil()
In [60]: Cl[r,c]=A.tolil()[r,c]
In [61]: Cl.tocsr()

C[r,c]=A[r,c] 给出了效率警告，但我认为这更多是为了让人们在循环中进行这种赋值。

In [63]: %%timeit C=A*B
    ...: C[r,c]=A[r,c]
...
The slowest run took 7.32 times longer than the fastest....
1000 loops, best of 3: 334 µs per loop

In [64]: %%timeit C=A*B
    ...: Cl=C.tolil()
    ...: Cl[r,c]=A.tolil()[r,c]
    ...: Cl.tocsr()
    ...: 
100 loops, best of 3: 2.83 ms per loop

我的 A 很小，只有 (250,100)，但看起来往返 lil 并不节省时间，尽管有警告。

当 A 稀疏时，使用 A==0 进行掩码必然会出现问题

In [66]: Az=A==0
....SparseEfficiencyWarning...
In [67]: r1,c1=Az.nonzero()

与 A 的 nonzero r 相比，这个 r1 大得多 - 所有的行索引稀疏矩阵中的零点；除了 25 个非零值以外的所有内容。

In [70]: r.shape
Out[70]: (25,)

In [71]: r1.shape
Out[71]: (24975,)

如果我用 r1 索引 A 我会得到一个更大的数组。实际上，我按其中零的数量重复每一行

In [72]: A[r1,:]
Out[72]: 
<24975x100 sparse matrix of type '<class 'numpy.float64'>'
    with 2473 stored elements in Compressed Sparse Row format>

In [73]: A
Out[73]: 
<250x100 sparse matrix of type '<class 'numpy.float64'>'
    with 25 stored elements in Compressed Sparse Row format>

我将非零元素的形状和数量增加了大约 100(列数)。

定义 foo，并复制 Divakar 的测试:

def foo(A,B):
    r,c = A.nonzero()
    C = A*B
    C[r,c] = A[r,c]
    return C

In [83]: timeit naive(A,B)
100 loops, best of 3: 2.53 ms per loop

In [84]: timeit proposed(A,B)
/...
  SparseEfficiencyWarning)
100 loops, best of 3: 4.48 ms per loop

In [85]: timeit foo(A,B)
...
  SparseEfficiencyWarning)
100 loops, best of 3: 2.13 ms per loop

所以我的版本有适度的速度改进。正如 Divakar 发现的那样，改变稀疏度会改变相对优势。我希望尺寸也会改变它们。

A.nonzero 使用 coo 格式的事实表明，使用该格式构建新数组可能是可行的。许多稀疏代码通过coo值构建一个新矩阵。

In [97]: Co=C.tocoo()    
In [98]: Ao=A.tocoo()

In [99]: r=np.concatenate((Co.row,Ao.row))
In [100]: c=np.concatenate((Co.col,Ao.col))
In [101]: d=np.concatenate((Co.data,Ao.data))

In [102]: r.shape
Out[102]: (79,)

In [103]: C1=sparse.csr_matrix((d,(r,c)),shape=A.shape)

In [104]: C1
Out[104]: 
<250x100 sparse matrix of type '<class 'numpy.float64'>'
    with 78 stored elements in Compressed Sparse Row format>

我认为，这个 C1 与通过其他方式构造的 C 具有相同的非零元素。但我认为一个值是不同的，因为 r 更长。在此特定示例中，C 和 A 共享一个非零元素，并且 coo 输入样式将这些元素相加，我们希望使用 A 值覆盖所有内容。

如果您可以容忍这种差异，这是一种更快的方法(至少对于这个测试用例而言):

def bar(A,B):
    C=A*B
    Co=C.tocoo()
    Ao=A.tocoo()
    r=np.concatenate((Co.row,Ao.row))
    c=np.concatenate((Co.col,Ao.col))
    d=np.concatenate((Co.data,Ao.data))
    return sparse.csr_matrix((d,(r,c)),shape=A.shape)

In [107]: timeit bar(A,B)
1000 loops, best of 3: 1.03 ms per loop

关于python - 两个仅影响零值的稀疏矩阵的点积，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/38419271/