在numpy中获得矩阵<矩阵结果的最快方法是什么？答案

【问题标题】：What is the fastest way to get the result of matrix < matrix in numpy?在numpy中获得矩阵<矩阵结果的最快方法是什么？
【发布时间】：2018-06-20 07:52:33
【问题描述】：

假设我有一个维度为 (M, A) 的矩阵 M_1 和一个维度为 (M, B) 的矩阵 M_2。 M_1 < M_2 的结果应该是一个维度为 (M, B, A) 的矩阵，其中 M1 中的每一行与 M_2 的相应行的每个元素进行比较，并给出一个布尔向量（或 1， 0-vector) 用于每次比较。

例如，如果我有一个矩阵

M1 = [[1,2,3]
      [3,4,5]]

M2 = [[1,2],
      [3,4]]

result should be [[[False, False, False],
                   [True, False, False]],
                  [[False, False, False], 
                   [True, False, False]]]

目前，我正在使用 for 循环，当我必须多次重复此操作（需要几个月）时，它的速度非常慢。希望有一种矢量化的方式来做到这一点。如果没有，我还能做什么？

我正在查看 M_1 是 (500, 3000000) 和 M_2 是 (500, 500) 并重复了大约 10000 次。

【问题讨论】：

@Divakar 抱歉已编辑
由于这是一个完全受内存限制的问题（在 numexpr 中，如@Divakar 所示，或者如果使用 Numba 或 Cython），最好知道之后你将如何处理这个布尔数组。
我不相信你可能真的需要这个。
另外一个 (500, 500, 3000000) 布尔矩阵几乎是 700 GiB。按位表示将其减少到 90 GiB 以下，但这意味着您不能使用基本的 Numpy 矩阵。
@Veedrac 然后让我向您解释一下设置。问题是在一个包含 300 万张人脸的数据库中测试人脸识别算法，平均约有 500 张人脸属于一个人/一个身份。测试是：对于每个身份，我需要计算属于同一身份（类内）的面孔之间的成对距离，并计算该身份的面孔与其他身份的面孔（类间）之间的成对距离。然后，理想情况下，我们应该有 max_intraclass_distance

标签： python performance numpy vectorization

【解决方案1】：

对于 NumPy 数组，使用 None/np.newaxis 扩展 dims 以便第一个轴对齐，而第二个轴是 spread 以便以元素方式比较它们。最后利用broadcasting 进行比较以获得矢量化解决方案 -

M1[:,None,:] < M2[:,:,None]

示例运行 -

In [19]: M1
Out[19]: 
array([[1, 2, 3],
       [3, 4, 5]])

In [20]: M2
Out[20]: 
array([[1, 2],
       [3, 4]])

In [21]: M1[:,None,:] < M2[:,:,None]
Out[21]: 
array([[[False, False, False],
        [ True, False, False]],

       [[False, False, False],
        [ True, False, False]]])

对于作为输入的列表，使用numpy.expand_dims，然后比较 -

In [42]: M1 = [[1,2,3],
    ...:       [3,4,5]]
    ...: 
    ...: M2 = [[1,2],
    ...:       [3,4]]

In [43]: np.expand_dims(M1, axis=1) < np.expand_dims(M2, axis=2)
Out[43]: 
array([[[False, False, False],
        [ True, False, False]],

       [[False, False, False],
        [ True, False, False]]])

进一步提升

利用 multi-core with numexpr module 处理大数据进一步提升 -

In [44]: import numexpr as ne

In [52]: M1 = np.random.randint(0,9,(500, 30000))

In [53]: M2 = np.random.randint(0,9,(500, 500))

In [55]: %timeit M1[:,None,:] < M2[:,:,None]
1 loop, best of 3: 3.32 s per loop

In [56]: %timeit ne.evaluate('M1e<M2e',{'M1e':M1[:,None,:],'M2e':M2[:,:,None]})
1 loop, best of 3: 1.53 s per loop

【讨论】：