给定两个 2D numpy 数组 A 和 B，如何有效地将采用两个 1D 数组的函数应用于 A 和 B 行的每个组合？

Posted 2023-02-14

【中文标题】给定两个 2D numpy 数组 A 和 B，如何有效地将采用两个 1D 数组的函数应用于 A 和 B 行的每个组合？

【英文标题】：Given two 2D numpy arrays A and B, how to efficiently apply a function that takes two 1D arrays to each combination of rows of A and B?

【发布时间】：2022-01-20 16:27:24

【问题描述】：

要清楚，以下是我想要做的。问题是，我如何更改函数oper_AB()，而不是嵌套的for循环，而是利用numpy中的矢量化/广播并更快地到达ret_list？

def oper(a_1D, b_1D):
    return np.dot(a_1D, b_1D) / np.dot(b_1D, b_1D)

def oper_AB(A_2D, B_2D):
    ret_list = []
    for a_1D in A_2D:
        for b_1D in B_2D:
            ret_list.append(oper(a_1D, b_1D))
    return ret_list

【问题讨论】：

旁注：数量oper(a, b)是什么意思？ a 在 b 上的标量投影会不同 (a @ b / np.linalg.norm(b))。当除以b @ b 时，您正在除以范数平方。

这是向量投影中的系数。 en.wikipedia.org/wiki/Vector_projection

标量投影为dot(a, b) / norm(b)，具有几何意义（向量间夹角的余弦）。 a 到b 的矢量投影为dot(a,b) / dot(b,b) * b，具有明显的几何解释。您计算的系数本身没有意义（没有乘以b）。如果你以后要乘以b，那很好。否则，我不明白它的含义。仅供参考。

这其实是有含义的。它只是投影向量的归一化版本。这是一个分数，表明a_1D 与b_1D 在b_1D 的方向上有多相似。如果为1，则表示a_1D在b_1D上的投影等于b_1D。

【参考方案1】：

这应该可以。

result = (np.matmul(A_2D, B_2D.transpose())/np.sum(B_2D*B_2D,axis=1)).flatten()

但是由于缓存利用率，第二个实现会更快。

def oper_AB(A_2D, B_2D):
    b_squared = np.sum(B_2D*B_2D,axis=1).reshape([-1,1])
    b_normalized = B_2D/b_squared
    del b_squared
    returned_val = np.matmul(A_2D,b_normalized.transpose())
    return returned_val.flatten()

del 只是在 B_2D 分配的内存太大时才会出现，（或者只是我习惯于使用多个 GB 数组）

编辑：根据 A_1D - B_1D 的要求

def oper2_AB(A_2D, B_2D):
    output = np.zeros([A_2D.shape[0]*B_2D.shape[0],A_2D.shape[1]],dtype=A_2D.dtype)
    for i in range(len(A_2D)):
        output[i*len(B_2D):(i+1)*len(B_2D)] = A_2D[i]-B_2D
    return output

【讨论】：

谢谢，这很快！当oper(a_1D,b_1D) 返回a_1D-b_1D 而不是np.dot(a_1D,b_1D)/np.dot(b_1D,b_1D) 时，我该如何做类似的事情？我们可以安全地假设向量 a_1D 和 b_1D 的大小是相同的。

刚刚用你的额外问题编辑了我的答案，如果解决方案有效，请务必投票。

好的，谢谢，我预计没有循环就没有办法做减法，看起来像这样。

有一种使用 repmat 的方法，但是它会占用更多的内存并且比循环版本执行得慢，所以我不写它。 （它让我的电脑崩溃了）

【参考方案2】：

严格解决这个问题（保留我怀疑 OP 想要规范，而不是规范平方，如下面的除数）：

r = a @ b.T / np.linalg.norm(b, axis=1)**2

例子：

np.random.seed(0)
a = np.random.randint(0, 10, size=(2,2))
b = np.random.randint(0, 10, size=(2,2))

然后：

>>> a
array([[5, 0],
       [3, 3]])

>>> b
array([[7, 9],
       [3, 5]])

>>> oper_AB(a, b)
[0.2692307692307692,
 0.4411764705882353,
 0.36923076923076925,
 0.7058823529411765]

>>> a @ b.T / np.linalg.norm(b, axis=1)**2
array([[0.26923077, 0.44117647],
       [0.36923077, 0.70588235]])

>>> np.ravel(a @ b.T / np.linalg.norm(b, axis=1)**2)
array([0.26923077, 0.44117647, 0.36923077, 0.70588235])

速度：

n, m = 1000, 100
a = np.random.uniform(size=(n, m))
b = np.random.uniform(size=(n, m))

orig = %timeit -o oper_AB(a, b)
# 2.73 s ± 11 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

new = %timeit -o np.ravel(a @ b.T / np.linalg.norm(b, axis=1)**2)
# 2.22 ms ± 33.3 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

orig.average / new.average
# 1228.78 (speedup)

我们的解决方案比原来的解决方案快 1200 倍。

正确性：

>>> np.allclose(np.ravel(a @ b.T / np.linalg.norm(b, axis=1)**2), oper_AB(a, b))
True

大型阵列的速度，与@Ahmed AEK 的解决方案相比：

n, m = 2000, 2000
a = np.random.uniform(size=(n, m))
b = np.random.uniform(size=(n, m))

new = %timeit -o np.ravel(a @ b.T / np.linalg.norm(b, axis=1)**2)
# 86.5 ms ± 484 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
other = %timeit -o AEK(a, b)  # Ahmed AEK's answer
# 102 ms ± 379 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

我们的解决方案速度提高了 15% :-)

【讨论】：

两种实现在我的机器上运行 40 次时实际上具有相同的平均运行时间，我认为减少 15% 的时间是靠运气或 CPU 上的一些硬件加速，你也不应该 永远永远当你的意思是 a*a 时使用 a**2 ，因为这个幂比向量乘法慢很多，它实际上比线性缩放慢。