此功能的完美 numpy 实现答案

【问题标题】：Perfect numpy implementation for this function此功能的完美 numpy 实现
【发布时间】：2020-05-26 12:49:46
【问题描述】：

这似乎是一个更直接的问题。最后我会概括一下。

我正在尝试在 numpy 中使用此功能。我已经成功地使用了嵌套的 for 循环，但我想不出一个 numpy 的方法来做到这一点。

我的实现方式：

bs = 10 # batch_size
nb = 8 # number of bounding boxes
nc = 15 # number of classes

bbox = np.random.random(size=(bs, nb, 4)) # model output bounding boxes

p = np.random.random(size=(bs, nb, nc)) # model output probability
p = softmax(p, axis=-1)

s_rand = np.random.random(size=(nc, nc))
s = (s_rand + s_rand.T)/2 # similarity matrix

pp = np.random.random(size=(bs, nb, nc)) # proposed probability
pp = softmax(pp, axis=-1)

first_term = 0
for b in range(nb):
    for b_1 in range(nb):
        if b_1 == b:
            continue
        for l in range(nc):
            for l_1 in range(nc):
                first_term += (s[l, l_1] * (pp[:, b, l] - pp[:, b_1, l_1])**2)
second_term = 0
for b in range(nb):
    for l in range(nc):
        second_term += (np.linalg.norm(s[l, :], ord=1) * (pp[:, b, l] - p[:, b, l])**2)
second_term *= nb

epsilon = 0.5
output = ((1 - epsilon) * first_term) + (epsilon * second_term)

为了完成任务，我已努力删除循环并改用np.tile 和np.repeat。但想不出可行的办法。

我曾尝试在 google 上搜索类似这样的练习，这些练习可以帮助我在 numpy 中学习此类转换，但没有成功。

【问题讨论】：

你有一些输入和一些预期的输出数据吗？
不，因为这是我的一个模型的损失函数。但我可以分享 numpy 数组的形状。将更新问题。
只是编一些数字？
是的..它有效。已经测试过了。所以，我只需要一个更好的选择。顺便说一句，更新了问题。
我猜内部循环可以替换为 .dot 操作（产品总和）

标签： python numpy machine-learning loss-function

【解决方案1】：

P_hat.shape 是 (B,L)，S.shape 是 (L,L)，P.shape 是 (B,L)。

array_before_sum = S[None,:,None,:]*(P_hat[:,:,None,None]- P_hat[None,None,:,:])**2
array_after_sum = array_before_sum.sum(axis=(1,3))
array_sum_again = (array_after_sum*(1-np.ones((B,B)))).sum()
first_term = (1-epsilon)*array_sum_again

second_term = epsilon*(B*np.abs(S).sum(axis=1)[None,:]*(P_hat - P)**2).sum()

【讨论】：

还没有检查过这个，但它似乎是正确的。交叉检查后将在此处更新。但是一个小问题。你是怎么做的？我的意思是，我想不出这样的方法。你在哪里练习/学习？请告诉我..我也想学习。
已验证。输出不同。因此，不是正确的。你能解释一下你的实现吗？ @L.铱
我只是以我的方式实现方程，所以我没有检查输出是否相同。对我来说，方程中的下标和 numpy 数组的索引之间的关系对于理解这个问题非常重要。熟悉广播规则也很有用。此外，数组的形状也很有用。
在我看来，从第一项中的四次求和我们知道我们将一个四维数组转换为一个数字。所以，我使用广播规则将 (B, L,1,1) 形状 P_hat 和 (1,1,B',L') 形状添加到 (B,L,B',L')四维数组和幂它的每个元素，B=B'，L=L'。然后，将 S 的形状从 (L, L') 转换为 (1, L, 1,L') 以便将其广播为 (B,L,B',L')。然后，我首先对 L 和 L' 维度求和，因为这个求和没有约束，并得到一个 (B,B') 数组。
继续。我发现在约束 b≠b' 下，求和不包括对角线中的元素。因此，我利用单位矩阵来解决它。这些只是我大脑中的过程。可能有问题。仅供参考。

【解决方案2】：

我认为你可以同时使用einsum

first_term = np.einsum('km, ijklm -> i', s, (pp[..., None, None] - pp[:, None, None, ...])**2 )
second_term = np.einsum('k, ijk -> i', np.linalg.norm(s, axis = 1), (pp - p)**2 )

现在有一个问题：如果 nb 和 nc 变大，first_term 中的 ijklm 张量将会变大。您可能应该分配它以获得 3 个较小的张量：

first_term =    np.einsum('km, ijk, ijk -> i', s, pp, pp)  +\
                np.einsum('km, ilm, ilm -> i', s, pp, pp)  -\
            2 * np.einsum('km, ijk, ilm -> i', s, pp, pp)

这利用了(a-b)**2 = a**2 + b**2 - 2ab 的优势，让您可以将问题分成三个部分，每个部分都可以通过点积一步完成

【讨论】：

我会试试你的方法。但是你能解释一下 einsum 是如何工作的吗？另外，这将如何给我带来速度优势？
我试过你的方法。您提供的输出与我通过相关代码获得的输出不同。我认为您缺少条件b != b_1。你能核对一下吗？
今天刚刚浪费了我所有的空闲时间，所以我无法真正调试。也不知道如何实现这种不平等。我会留下来，希望能激发这里的其他einsum 大师之一。
没问题.. 感谢您的尝试。

【解决方案3】：

最大限度优化的代码：（删除前两个循环的灵感来自 L.Iridium 的回答）

squared_diff = (pp[:, :, None, :, None] - pp[:, None, :, None, :]) ** 2
weighted_diff = s * squared_diff
b_eq_b_1_removed = b.sum(axis=(3,4)) * (1 - np.eye(nb))
first_term = b_eq_b_1_removed.sum(axis=(1,2))

normalized_s = np.linalg.norm(s, ord=1, axis=1)
squared_diff = (pp - p)**2
second_term = nb * (normalized_s * squared_diff).sum(axis=(1,2))

loss = ((1 - epsilon) * first_term) + (epsilon * second_term)

时间轨迹： 512 µs ± 13 µs per loop

有问题的代码的Timeit跟踪： 62.5 ms ± 197 µs per loop

这是一个巨大的进步。

【讨论】：