Sparse Scipy/Numpy：一种实现两两分钟求和运算的有效方法答案

【问题标题】：Sparse Scipy/Numpy: an efficient way to implement sum of pairwise mins operationSparse Scipy/Numpy：一种实现两两分钟求和运算的有效方法
【发布时间】：2020-02-27 20:36:47
【问题描述】：

计算向量之间的成对最小值之和在自然语言处理 (NLP) 中非常流行，并用于计算相交直方图内核 [1]。然而，在 NLP 中，我们经常处理稀疏矩阵。

这是一种使用慢速 for 循环来计算此操作的低效方法：

import numpy as np
from scipy.sparse import csr_matrix

# Initialize sparse matrices
A = csr_matrix(np.clip(np.random.randn(100, 64) - 1, 0, np.inf))
B = csr_matrix(np.clip(np.random.randn(64, 100) - 1, 0, np.inf))

# For each row, col vector i,j in A and B respectively
G = np.zeros((100, 100))
for i in range(A.shape[0]):
    for j in range(B.shape[1]):
        G[i, j] = A[i].minimum(B[:,j]).sum()

有没有办法在没有 for 循环的情况下做到这一点？如果可以编译 for 循环，例如在 numba 中使用 jit，我不介意。

这里给出了一个快速密集的版本：Numpy: an efficient way to implement sum of pairwise mins operation

谢谢。

[1]http://blog.datadive.net/histogram-intersection-for-change-detection/

【问题讨论】：

非零值是否总是正数？
所以，它们并不是很稀疏。为什么要使用稀疏表示？是内存效率困扰你吗？
@Divakar 这只是为了说明，我使用的矩阵实际上非常稀疏。
@hpaulj，非零值都是正数！
那么，您的意思是它们是 NumPy 数组，但非常“稀疏”，就像许多零一样？

标签： python numpy scipy sparse-matrix

【解决方案1】：

这是一个应该是高效的实现，尽可能地利用稀疏性。有一个循环，但只有一个暗处，所以应该不会太糟糕。

import numpy as np
from scipy.sparse import csr_matrix, csc_matrix

M, N, K = 640, 100, 650

B1 = csr_matrix(np.clip(np.random.randn(N, K) - 1, 0, np.inf))
B2 = csr_matrix(np.clip(np.random.randn(N, K) - 1, 0, np.inf))
B = B1-B2
A1 = csc_matrix(np.clip(np.random.randn(M, N) - 1, 0, np.inf))
A2 = csc_matrix(np.clip(np.random.randn(M, N) - 1, 0, np.inf))
A = A1-A2

result = np.zeros((M, K))
for j in range(N):
    ia = A.indices[A.indptr[j] : A.indptr[j+1]]
    ib = B.indices[B.indptr[j] : B.indptr[j+1]]
    IA, IB = np.ix_(ia, ib)
    da = A.data[A.indptr[j] : A.indptr[j+1]]
    db = B.data[B.indptr[j] : B.indptr[j+1]]
    # both nonzero
    result[IA, IB] += np.minimum.outer(da, db)
    # one negative ...
    am = da<0
    iam, dam = ia[am], da[am]
    bm = db<0
    ibm, dbm = ib[bm], db[bm]
    # ... the other zero
    za = np.ones((M,), dtype=bool)
    za[ia] = False
    zb = np.ones((K,), dtype=bool)
    zb[ib] = False
    IA, IB = np.ix_(iam, zb)
    result[IA, IB] += dam[:, None]
    IA, IB = np.ix_(za, ibm)
    result[IA, IB] += dbm
# compare with dense method
print(np.allclose(result, np.minimum(A.A[..., None], B.A).sum(axis=1)))

打印

True

【讨论】：

看起来不错，打算试试看它是否与numba一起编译
@Curious 很高兴听到它。顺便提一句。如果您的矩阵没有负条目，则不需要从# one negative ... 到循环末尾的代码。你可以简单地省略它。事实上，你应该这样做，因为这会让事情变得更快。

【解决方案2】：

好吧，至少在 SciPy 的最新版本中，有一个函数 scipy.sparse.csr_matrix.minimum Link to documentation 在元素最小值方面相当于 numpy.minimum。但是，我不知道计算效率如何。

【讨论】：