Numpy数组比列表慢得多

Question

给定两个矩阵X1（N，3136）和X2（M，3136）（其中每行中的每个元素都是二进制数），我试图计算汉明距离，以便将X1中的每个元素与X2的所有行进行比较，结果矩阵为（N，M）。

我编写了两个函数来实现此功能（第一个使用numpy帮助，另一个不使用numpy）：

def hamming_distance(X, X_train):
    array = np.array([np.sum(np.logical_xor(x, X_train), axis=1) for x in X])

    return array



def hamming_distance2(X, X_train):
    a = len(X[:,0])
    b = len(X_train[:,0])

    hamming_distance = np.zeros(shape=(a, b))

    for i in range(0, a):
        for j in range(0, b):
            hamming_distance[i,j] = np.count_nonzero(X[i,:] != X_train[j,:])

    return hamming_distance

我的问题是，使用两个for循环的下面的函数比使用upper函数的上面的函数慢得多。有没有可能改进第一个函数，以便我只使用一个循环？

PS. 对不起，我的英语不是我的母语，尽管我已经尽力了！

- Derik Daro

---

从中学习（(u != v).mean()）。 - sascha

“(u != v).mean()”不会比当前版本更有效率，因为循环仍然是必要的。N和M有多大？也许可以通过广播计算所有成对比较，并仅切片需要的部分。但这可能会很容易地填满您的内存。 - ayhan

这只是说明numpy不是魔法灰尘；只有在适当使用时才能提高性能。而你的第二个也使用了numpy。 - Andras Deak -- Слава Україні

我的失败尝试是使用点积，因为这两个操作非常相似。给定两个二进制矩阵，如果你将其中一个与(1-另一个)相乘，它会给出不同条目的数量，其中第一个是1，而另一个是0。所以(1-arr1).dot(arr2.T) + arr1.dot(1-arr2.T)实际上给出了汉明距离。但是，它比循环更快吗？不幸的是，不是。 :) - ayhan

为了提取最佳性能，我们需要利用带有浮点值的BLAS。请发布以float32转换数组为基础的点积答案帖子。它的胜利幅度很大。 - Divakar

Answer 1

Numpy只有在你使用向量化处理时才能使你的代码更快。在您的情况下，您可以利用数组广播来向量化您的问题：比较您的两个数组并创建一个形状为(N,M,K)的辅助数组，您可以沿其第三维进行求和。

hamming_distance = (X[:,None,:] != X_train).sum(axis=-1)

我们向第一个数组中注入一个单例维度，使其形状为(N,1,K)，第二个数组隐式兼容形状(1,M,K)，因此可以执行操作。

在评论中@ayhan指出，这将为大型M和N创建一个巨大的辅助数组，这是非常正确的。这就是向量化的代价：你通过内存来获得CPU时间。如果你有足够的内存让上述方法运行，它会非常快。如果没有，你必须减少向量化的范围，并在M或N（或两者）中循环（这将成为您当前的方法）。但这与numpy本身无关，这涉及在可用资源和性能之间取得平衡。

- Andras Deak -- Слава Україні

---

这将生成大约470亿个条目。 - ayhan

@ayhan 我猜那是正确的 :D 没有仔细阅读评论。嗯，这就是向量化所需要付出的代价。 - Andras Deak -- Слава Україні

@ayhan 现在我仔细阅读了评论，我刚刚注意到你基本上提出了相同的解决方案。你想加入自己的答案吗？如果是这样，我很乐意删除我的答案。 - Andras Deak -- Слава Україні

没错，当我看到那些数字的时候，我就退缩了。:) 别担心，我正在努力寻找另一个解决方案。 - ayhan

Answer 2

您正在做的事情非常类似于点积。考虑这两个二进制数组：

1 0 1 0 1 1 0 0
0 0 1 1 0 1 0 1

我们正在努力寻找不同对数。如果您直接取点积，它将为您提供（1，1）对的数量。然而，如果你否定其中一个，它会计算不同的数量。例如，a1.dot(1-a2) 计算 (1, 0) 对。由于我们还需要 (0, 1) 对的数量，因此我们将 a2.dot(1-a1) 加入其中。点积的好处是它非常快速。但是，您需要先将数组转换为浮点数，正如Divakar指出的那样。

以下是演示：

prng = np.random.RandomState(0)
arr1 = prng.binomial(1, 0.3, (1000, 3136))
arr2 = prng.binomial(1, 0.3, (2000, 3136))
res1 = hamming_distance2(arr1, arr2)
arr1 = arr1.astype('float32'); arr2 = arr2.astype('float32')
res2 = (1-arr1).dot(arr2.T) + arr1.dot(1-arr2.T)

np.allclose(res1, res2)
Out: True

还有时间安排：

%timeit hamming_distance(arr1, arr2)
1 loop, best of 3: 13.9 s per loop

%timeit hamming_distance2(arr1, arr2)
1 loop, best of 3: 5.01 s per loop

%timeit (1-arr1).dot(arr2.T) + arr1.dot(1-arr2.T)
10 loops, best of 3: 93.1 ms per loop

- ayhan