两个稀疏矩阵的点积仅影响零值

Question

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

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太远.我的想法是将dot产品保留为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会有所帮助.

Answer 1

代码的可能更干净、更快的版本naive：

\n\n

In [57]: r,c=A.nonzero()    # this uses A.tocoo()\n\nIn [58]: C=A*B\nIn [59]: Cl=C.tolil()\nIn [60]: Cl[r,c]=A.tolil()[r,c]\nIn [61]: Cl.tocsr()\n

\n\n

C[r,c]=A[r,c]给出了效率警告，但我认为这更多的是为了让人们在循环中完成这种任务。

\n\n

In [63]: %%timeit C=A*B\n    ...: C[r,c]=A[r,c]\n...\nThe slowest run took 7.32 times longer than the fastest....\n1000 loops, best of 3: 334 \xc2\xb5s per loop\n\nIn [64]: %%timeit C=A*B\n    ...: Cl=C.tolil()\n    ...: Cl[r,c]=A.tolil()[r,c]\n    ...: Cl.tocsr()\n    ...: \n100 loops, best of 3: 2.83 ms per loop\n

\n\n

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

\n\n

当稀疏时，掩蔽A==0必然会出现问题A

\n\n

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

\n\n

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

\n\n

In [70]: r.shape\nOut[70]: (25,)\n\nIn [71]: r1.shape\nOut[71]: (24975,)\n

\n\n

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

\n\n

In [72]: A[r1,:]\nOut[72]: \n<24975x100 sparse matrix of type '<class 'numpy.float64'>'\n    with 2473 stored elements in Compressed Sparse Row format>\n\nIn [73]: A\nOut[73]: \n<250x100 sparse matrix of type '<class 'numpy.float64'>'\n    with 25 stored elements in Compressed Sparse Row format>\n

\n\n

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

\n\n

定义foo并复制 Divakar 的测试：

\n\n

def foo(A,B):\n    r,c = A.nonzero()\n    C = A*B\n    C[r,c] = A[r,c]\n    return C\n\nIn [83]: timeit naive(A,B)\n100 loops, best of 3: 2.53 ms per loop\n\nIn [84]: timeit proposed(A,B)\n/...\n  SparseEfficiencyWarning)\n100 loops, best of 3: 4.48 ms per loop\n\nIn [85]: timeit foo(A,B)\n...\n  SparseEfficiencyWarning)\n100 loops, best of 3: 2.13 ms per loop\n

\n\n

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

\n\n

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

\n\n

In [97]: Co=C.tocoo()    \nIn [98]: Ao=A.tocoo()\n\nIn [99]: r=np.concatenate((Co.row,Ao.row))\nIn [100]: c=np.concatenate((Co.col,Ao.col))\nIn [101]: d=np.concatenate((Co.data,Ao.data))\n\nIn [102]: r.shape\nOut[102]: (79,)\n\nIn [103]: C1=sparse.csr_matrix((d,(r,c)),shape=A.shape)\n\nIn [104]: C1\nOut[104]: \n<250x100 sparse matrix of type '<class 'numpy.float64'>'\n    with 78 stored elements in Compressed Sparse Row format>\n

\n\n

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

\n\n

如果您可以容忍这种差异，这是一种更快的方法（至少对于这个测试用例来说）：

\n\n

def bar(A,B):\n    C=A*B\n    Co=C.tocoo()\n    Ao=A.tocoo()\n    r=np.concatenate((Co.row,Ao.row))\n    c=np.concatenate((Co.col,Ao.col))\n    d=np.concatenate((Co.data,Ao.data))\n    return sparse.csr_matrix((d,(r,c)),shape=A.shape)\n\nIn [107]: timeit bar(A,B)\n1000 loops, best of 3: 1.03 ms per loop\n

\n