计算略有不同的矩阵乘法答案

【问题标题】：Computing a slightly different matrix multiplication计算略有不同的矩阵乘法
【发布时间】：2019-06-05 10:48:17
【问题描述】：

我正在尝试找到计算两组向量之间最小元素乘积的最佳方法。通常的矩阵乘法C=A@B 将Cij 计算为向量Ai 和B^Tj 的元素的成对乘积之和。我想执行成对产品中的最小值。我找不到一个有效的方法来在两个矩阵之间使用 numpy.

实现此目的的一种方法是生成A 和B（在求和之前）之间的成对乘积的 3D 矩阵，然后在第三维上取最小值。但这会导致巨大的内存占用（实际上我不知道该怎么做）。

你知道我如何实现这个操作吗？

示例：

A = [[1,1],[1,1]]
B = [[0,2],[2,1]]

矩阵matmul：

C = [[1*0+1*2,1*2+1*1][1*0+1*2,1*2+1*1]] = [[2,3],[2,3]]

最小matmul：

C = [[min(1*0,1*2),min(1*2,1*1)][min(1*0,1*2),min(1*2,1*1)]] = [[0,1],[0,1]]

【问题讨论】：

标签： python numpy scipy sparse-matrix

【解决方案1】：

在将A 扩展为3D 后使用broadcasting -

A = np.asarray(A)
B = np.asarray(B)
C_out = np.min(A[:,None]*B,axis=2)

如果您关心内存占用，请使用numexpr module 以提高效率 -

import numexpr as ne

C_out = ne.evaluate('min(A3D*B,2)',{'A3D':A[:,None]})

大型数组的计时 -

In [12]: A = np.random.rand(200,200)

In [13]: B = np.random.rand(200,200)

In [14]: %timeit np.min(A[:,None]*B,axis=2)
34.4 ms ± 614 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [15]: %timeit ne.evaluate('min(A3D*B,2)',{'A3D':A[:,None]})
29.3 ms ± 316 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [16]: A = np.random.rand(300,300)

In [17]: B = np.random.rand(300,300)

In [18]: %timeit np.min(A[:,None]*B,axis=2)
113 ms ± 2.27 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [19]: %timeit ne.evaluate('min(A3D*B,2)',{'A3D':A[:,None]})
102 ms ± 691 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

所以，numexpr 有一些改进，但可能没有我预期的那么好。

【讨论】：

很好，您能详细说明一下 numexpr 模块吗？该文档不是很有帮助。这是否意味着将通过计算每对向量的最小值来评估表达式，而不存储实际的成对乘积？
@debzsud 我不知道它的确切内部实现，但从文档中它声明它使用多核以及从以前的解决方案的基准测试中，我已经看到它在大型阵列上的性能非常好，而且最肯定应该是内存效率。
@Divakar 在这个例子中你观察到多核的使用吗？
@max9111 添加了时间和结论。内存分析器可以进一步澄清事情。
你知道如何应用广播版本（第一个版本），A 和 B 是 scipy.sparse.csr_matrix 实例吗？

【解决方案2】：

Numba 也是一种选择

我对不是特别好的 Numexpr Timings 有点惊讶，所以我尝试了 Numba 版本。对于大型阵列，这可以进一步优化。（可以应用与 dgemm 相同的原理）

import numpy as np
import numba as nb
import numexpr as ne

@nb.njit(fastmath=True,parallel=True)
def min_pairwise_prod(A,B):
    assert A.shape[1]==B.shape[1]
    res=np.empty((A.shape[0],B.shape[0]))

    for i in nb.prange(A.shape[0]):
        for j in range(B.shape[0]):
            min_prod=A[i,0]*B[j,0]
            for k in range(B.shape[1]):
                prod=A[i,k]*B[j,k]
                if prod<min_prod:
                    min_prod=prod
            res[i,j]=min_prod
    return res

时间安排

A=np.random.rand(300,300)
B=np.random.rand(300,300)

%timeit res_1=min_pairwise_prod(A,B) #parallel=True
5.56 ms ± 1.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
%timeit res_1=min_pairwise_prod(A,B) #parallel=False
26 ms ± 163 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
%timeit res_2 = ne.evaluate('min(A3D*B,2)',{'A3D':A[:,None]})
87.7 ms ± 265 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
%timeit res_3=np.min(A[:,None]*B,axis=2)
110 ms ± 214 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

A=np.random.rand(1000,300)
B=np.random.rand(1000,300)
%timeit res_1=min_pairwise_prod(A,B) #parallel=True
50.6 ms ± 401 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)
%timeit res_1=min_pairwise_prod(A,B) #parallel=False
296 ms ± 5.02 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
%timeit res_2 = ne.evaluate('min(A3D*B,2)',{'A3D':A[:,None]})
992 ms ± 7.59 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
%timeit res_3=np.min(A[:,None]*B,axis=2)
1.27 s ± 15.4 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

【讨论】：