python - bool 矩阵计算的最快方法

Question

我有一个包含 1.5E6 行和 20E3 列的 bool 矩阵，类似于以下示例:

M = [[ True,  True, False,  True, ...],
     [False,  True,  True,  True, ...],
     [False, False, False, False, ...],
     [False,  True, False, False, ...],
     ...
     [ True,  True, False, False, ...]
     ]

另外，我还有另一个矩阵 N ( 1.5E6行， 1列):

 N = [[ True],
      [False],
      [ True],
      [ True],
      ...
      [ True]
      ]

我需要做的是通过 M运算符组合的矩阵 AND(1＆1、1、2、1＆3、1＆N，2＆1、2＆2等)中的每一列对，并计算结果和矩阵 N之间有多少重叠。

我的Python/Numpy代码如下所示:

for i in range(M.shape[1]):
  for j in range(M.shape[1]):
    result = M[:,i] & M[:,j] # Combine the columns with AND operator
    count = np.sum(result & N.ravel()) # Counts the True occurrences
    ... # Save the count for the i and j pair

问题是，要通过带有两个for循环的 20E3 x 20E3组合在计算上是昂贵的( 需要大约5-10天的时间才能计算出)。我尝试过的一个更好的选择是将每一列与整个矩阵M进行比较:

for i in range(M.shape[1]):
  result = M[:,i]*M.shape[1] & M # np.tile or np.repeat is used to horizontally repeat the column
  counts = np.sum(result & N*M.shape[1], axis=0)
  ... # Save the counts

这样可以将开销和计算时间减少到10％左右，但是 仍需要1天左右的时间来计算。

我的问题是:进行这些计算(基本上只是 AND和 SUM)的最快方法(也许不是Python？)是什么？

我当时在考虑低级语言，GPU处理，量子计算等。但是我对这些都不了解，因此对方向的任何建议都非常感谢!

其他想法:
当前正在考虑是否可以使用点积(如Davikar提出的那样)快速计算组合的三元组:

def compute(M, N):
    out = np.zeros((M.shape[1], M.shape[1], M.shape[1]), np.int32)
    for i in range(M.shape[1]):
        for j in range(M.shape[1]):
            for k in range(M.shape[1]):
                result = M[:, i] & M[:, j] & M[:, k]
                out[i, j, k] = np.sum(result & N.ravel())
    return out

Answer 1

只需使用 np.einsum 即可获取所有计数-

np.einsum('ij,ik,i->jk',M,M.astype(int),N.ravel())

随意使用带有 optimize的 np.einsum标志。此外，还可以随意使用不同的dtypes转换。

要利用GPU，我们可以使用 tensorflow包，该包也支持 einsum 。

使用np.dot的更快替代品:

(M&N).T.dot(M.astype(int))
(M&N).T.dot(M.astype(np.float32))

时间-

In [110]: np.random.seed(0)
     ...: M = np.random.rand(500,300)>0.5
     ...: N = np.random.rand(500,1)>0.5

In [111]: %timeit np.einsum('ij,ik,i->jk',M,M.astype(int),N.ravel())
     ...: %timeit (M&N).T.dot(M.astype(int))
     ...: %timeit (M&N).T.dot(M.astype(np.float32))
227 ms ± 191 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)
66.8 ms ± 198 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
3.26 ms ± 753 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

并通过两个 bool 数组的float32转换使它更进一步-

In [122]: %%timeit
     ...: p1 = (M&N).astype(np.float32)
     ...: p2 = M.astype(np.float32)
     ...: out = p1.T.dot(p2)
2.7 ms ± 34.3 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)