python - 使用 scipy.sparse.csc_matrix 替代 numpy 广播

Question

我的代码中有以下表达式:

a = (b / x[:, np.newaxis]).sum(axis=1)

其中 b 是形状为 (M, N) 的数组，x 是形状为 (M, )。现在，b 实际上是稀疏的，因此为了提高内存效率，我想用 scipy.sparse.csc_matrix 或 csr_matrix 代替。但是，以这种方式进行的广播并未实现(即使保证除法或乘法保持稀疏性)(x 的条目非零)，并引发 NotImplementedError .是否有一个我不知道的 sparse 函数可以满足我的要求？ (dot() 会沿着错误的轴求和。)

Answer 1

如果b 是CSC 格式，那么b.data 有b 和b 的非零条目。 indices 具有每个非零条目的行索引，因此您可以按以下方式进行除法:

b.data /= np.take(x, b.indices)

它比 Warren 的优雅解决方案更 hack，但在大多数设置中它也可能更快:

b = sps.rand(1000, 1000, density=0.01, format='csc')
x = np.random.rand(1000)

def row_divide_col_reduce(b, x):
    data = b.data.copy() / np.take(x, b.indices)
    ret = sps.csc_matrix((data, b.indices.copy(), b.indptr.copy()),
                         shape=b.shape)
    return ret.sum(axis=1)

def row_divide_col_reduce_bis(b, x):
    d = sps.spdiags(1.0/x, 0, len(x), len(x))
    return (d * b).sum(axis=1)

In [2]: %timeit row_divide_col_reduce(b, x)
1000 loops, best of 3: 210 us per loop

In [3]: %timeit row_divide_col_reduce_bis(b, x)
1000 loops, best of 3: 697 us per loop

In [4]: np.allclose(row_divide_col_reduce(b, x),
   ...:             row_divide_col_reduce_bis(b, x))
Out[4]: True

如果您就地进行除法，则可以将上述示例中的时间几乎减半，即:

def row_divide_col_reduce(b, x):
    b.data /= np.take(x, b.indices)
    return b.sum(axis=1)

In [2]: %timeit row_divide_col_reduce(b, x)
10000 loops, best of 3: 131 us per loop