Numpy广播以矢量形式执行欧氏距离

Question

我有2×4和3×4的矩阵。我想找出横排的欧几里德距离，在最后得到一个2×3矩阵。这是一个for循环的代码，它根据所有b行向量计算a中每行向量的欧几里德距离。在不使用for循环的情况下如何执行相同的操作？

 import numpy as np
a = np.array([[1,1,1,1],[2,2,2,2]])
b = np.array([[1,2,3,4],[1,1,1,1],[1,2,1,9]])
dists = np.zeros((2, 3))
for i in range(2):
      dists[i] = np.sqrt(np.sum(np.square(a[i] - b), axis=1))

Answer 1

我最近在深度学习方面也遇到了同样的问题（斯坦福cs231n，作业1），但是当我使用

 np.sqrt((np.square(a[:,np.newaxis]-b).sum(axis=2)))

有个错误

MemoryError

这意味着我的内存不足（事实上，这中间产生了一个500*5000*1024的数组，太大了！）

为了防止这种错误，我们可以使用一个公式来简化：

代码：

import numpy as np
aSumSquare = np.sum(np.square(a),axis=1);
bSumSquare = np.sum(np.square(b),axis=1);
mul = np.dot(a,b.T);
dists = np.sqrt(aSumSquare[:,np.newaxis]+bSumSquare-2*mul)

Answer 2

只需在正确的位置使用np.newaxis：

 np.sqrt((np.square(a[:,np.newaxis]-b).sum(axis=2)))

Answer 3

以下是原始输入变量：

A = np.array([[1,1,1,1],[2,2,2,2]])
B = np.array([[1,2,3,4],[1,1,1,1],[1,2,1,9]])
A
# array([[1, 1, 1, 1],
#        [2, 2, 2, 2]])
B
# array([[1, 2, 3, 4],
#        [1, 1, 1, 1],
#        [1, 2, 1, 9]])

A是2x4阵列。 B是一个3x4阵列。

我们想在一个完全矢量化的运算中计算欧几里德距离矩阵运算，其中dist[i,j]包含A中的第i个实例和B中的第j个实例之间的距离。因此在这个例子中dist是2x3。

距离

表面上可以用numpy写成

dist = np.sqrt(np.sum(np.square(A-B))) # DOES NOT WORK
# Traceback (most recent call last):
#   File "<stdin>", line 1, in <module>
# ValueError: operands could not be broadcast together with shapes (2,4) (3,4)

然而，如上所示，问题在于按元素进行的减法操作A-B涉及不兼容的数组大小，特别是第一维中的2和3。

A has dimensions 2 x 4
B has dimensions 3 x 4

为了进行元素相减，我们必须填充A或B以满足numpy的广播规则。我将选择用一个额外的维度来填充，这样它就变成2 x 1 x 4，这样数组的维度就可以排列起来进行广播。有关numpy广播的更多信息，请参见tutorial in the scipy manual和this tutorial中的最后一个示例。

您可以使用np.newaxis值或使用np.reshape命令执行填充。我在下面展示两个：

# First approach is to add the extra dimension to A with np.newaxis
A[:,np.newaxis,:] has dimensions 2 x 1 x 4
B has dimensions                     3 x 4

# Second approach is to reshape A with np.reshape
np.reshape(A, (2,1,4)) has dimensions 2 x 1 x 4
B has dimensions                          3 x 4

如您所见，使用任何一种方法都将允许维度对齐。我将对np.newaxis使用第一种方法。现在，这将创建一个A-B，它是一个2x3x4数组：

diff = A[:,np.newaxis,:] - B
# Alternative approach:
# diff = np.reshape(A, (2,1,4)) - B
diff.shape
# (2, 3, 4)

现在我们可以将这个差分表达式放入dist等式语句中，得到最终结果：

dist = np.sqrt(np.sum(np.square(A[:,np.newaxis,:] - B), axis=2))
dist
# array([[ 3.74165739,  0.        ,  8.06225775],
#        [ 2.44948974,  2.        ,  7.14142843]])

注意sum在axis=2之上，这意味着取2x3x4数组第三个轴上的和（其中轴id以0开头）。

如果数组很小，那么上面的命令就可以正常工作。但是，如果有大型数组，则可能会遇到内存问题。注意，在上面的例子中，numpy内部创建了一个2x3x4数组来执行广播。如果我们将A推广到有维度a x z，B推广到有维度b x z，那么numpy将在内部创建一个用于广播的a x b x z数组。

我们可以通过一些数学操作来避免创建这个中间数组。因为你把欧几里德距离计算成平方差之和，我们可以利用平方差之和可以重写的数学事实。

注意，中间项涉及元素的和乘法。这个乘法和被称为点积。因为A和B都是A矩阵，所以这个运算实际上是矩阵乘法。因此，我们可以将上述内容改写为：

然后我们可以编写以下numpy代码：

threeSums = np.sum(np.square(A)[:,np.newaxis,:], axis=2) - 2 * A.dot(B.T) + np.sum(np.square(B), axis=1)
dist = np.sqrt(threeSums)
dist
# array([[ 3.74165739,  0.        ,  8.06225775],
#        [ 2.44948974,  2.        ,  7.14142843]])

请注意，上面的答案与前面的实现完全相同。同样，这里的优势是我们不需要创建用于广播的中间2x3x4阵列。

为了完整性，让我们仔细检查一下threeSums中每个summand的维度是否允许广播。

np.sum(np.square(A)[:,np.newaxis,:], axis=2) has dimensions 2 x 1
2 * A.dot(B.T) has dimensions                               2 x 3
np.sum(np.square(B), axis=1) has dimensions                 1 x 3

因此，正如预期的那样，最终的dist数组的维数是2x3。

在this tutorial中还讨论了使用点积代替按元素乘法的和的问题。