不同形状阵列的numpy距离计算答案

【问题标题】：Numpy distance calculations of different shaped arrays不同形状阵列的numpy距离计算
【发布时间】：2018-12-17 09:51:10
【问题描述】：

不确定我的标题是否正确，但基本上我有一个参考坐标，格式为 (x,y,z)，还有一个大的坐标列表/数组也采用该格式。我需要得到每个之间的欧几里得距离，所以理论上使用 numpy 和 scipy 我应该能够执行以下操作：

import numpy, scipy.spatial.distance
a = numpy.array([1,1,1])
b = numpy.random.rand(20,3)

distances = scipy.spatial.distance.euclidean(b, a)

但我没有得到一个数组，而是得到一个错误：ValueError: Input vector should be 1-D.

不知道如何解决这个错误并获得我想要的东西而不必求助于循环等，这有点违背了使用 Numpy 的目的。

从长远来看，我想使用这些距离来计算真值掩码，以计算 bin 中的距离值。

我不确定我只是使用了错误的函数还是使用了错误的函数，我无法在文档中找到任何可以更好地工作的内容。

【问题讨论】：

查看相关帖子stackoverflow.com/q/27948363/2588210
并查看文档：docs.scipy.org/doc/scipy/reference/generated/…
您正在阅读哪些文档？ docs.scipy.org/doc/scipy-1.0.0/reference/spatial.distance.html 的文档似乎很清楚这个函数需要一维输入。
如果你可以使用编译器，这样的问题真的很容易实现。例如。优于 cdist stackoverflow.com/a/49573091/4045774 还要准确输入您的输入内容（列表与数组非常不同）。关于 bincounting/digitize stackoverflow.com/a/50496823/4045774 的一些示例
@Will 以clang为例。它将 C 代码转换为 LLVM-IR，在此步骤之后，LLVM 后端生成机器代码。例如，Numba 将 Python 代码转换为 LLVM-IR，LLVM-IR 将由 LLVM 后端编译为机器代码。 Cython 将 Python 代码转换为 C，机器代码由 C 编译器创建。 Numba 仅支持一小部分 Python 语言或非常小的模块子集，但它是提高关键代码路径性能的好工具。

标签： python arrays numpy scipy euclidean-distance

【解决方案1】：

scipy.spatial.distance.euclidean 的 documentation 声明，仅允许一维向量作为输入。因此，您必须像这样遍历数组：

distances = np.empty(b.shape[0])
for i in range(b.shape[0]):
    distances[i] = scipy.spatial.distance.euclidean(a, b[i])

如果你想有一个矢量化的实现，你需要编写你自己的函数。也许使用带有正确签名的np.vectorize 也可以，但实际上这也只是 for 循环的简写，因此与简单的 for 循环具有相同的性能。

正如我对 hannes wittingham 解决方案的评论中所述，我将发布一个专注于性能的单行：

distances = ((b - a)**2).sum(axis=1)**0.5

写出所有计算减少了单独函数调用的数量，从而减少了将中间结果分配给新数组的次数。因此，对于b.shape == (20, 3) 的阵列形状，它比使用 hannes wittingham 的解决方案快约 22%，对于阵列形状为 b.shape == (20000, 3):

a = np.array([1, 1, 1,])
b = np.random.rand(20, 3)
%timeit ((b - a)**2).sum(axis=1)**0.5
# 5.37 µs ± 140 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
%timeit euclidean_distances(a, b)
# 6.89 µs ± 345 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

b = np.random.rand(20000, 3)
%timeit ((b - a)**2).sum(axis=1)**0.5
# 588 µs ± 43.2 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)
%timeit euclidean_distances(a, b)
# 616 µs ± 36.3 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

但是您正在失去能够轻松更改为距离计算例程的灵活性。使用scipy.spatial.distance模块时，只需调用其他方法即可更改计算路由。

为了进一步提高计算性能，您可以为您的函数使用像 numba 这样的 jit（即时）编译器：

import numba as nb
@nb.njit
def euc(a, b):
    return ((b - a)**2).sum(axis=1)**0.5

这将小型阵列的计算时间减少了大约 70%，而大型阵列的计算时间减少了大约 60%。不幸的是，np.linalg.norm 的 axis 关键字尚不被 numba 支持。

【讨论】：

我非常喜欢这个，我明白你关于失去灵活性的观点。我发现如果我使用 numpy.linalg 我可以做我正在寻找的事情，就像distances = numpy.linalg(a - b, axis=1) 这让我得到了我所期望的，它肯定不如 scipy 方法灵活，但它是一条相当快的路线。所以我认为这将是我使用一个或另一个的情况。单独贴出来供大家参考。
我猜你的意思是numpy.linalg.norm？这确实是进行计算的好方法。当我的代码不是性能关键时，我也会使用它（嗯，在 Python 中它通常不应该，但节省一些时间总是很好:)）。不知道为什么我忘了在我的回答中提到它。 np.linalg.norm 应该比 ((b - a)**2).sum(axis=1)**0.5 慢 5%-30%，具体取决于数组大小。非常感谢您标记我的答案，我很高兴能够提供帮助。
@scotty 以向量化的方式编写 Numba 函数并不值得推荐。有时它的工作速度与简单循环一样快（如果 Numba 可以将其翻译成这个），但很多时候它更慢。 @nb.njit(fastmath=True) 也可能是有益的（通常是 SIMD 矢量化所必需的）
没错，但由于这个问题不是关于 numba 或高度优化代码性能的问题，我只想发布一个简单而简短的示例来说明编译的含义。在这种情况下，与额外的努力相比，性能提升可以忽略不计。
@Scotty1- 这是一个很好的观点，30% 是计算时间的一个很大变化，我必须试一试，看看它是否能让我更快地运行，因为我的数据集是相当大。

【解决方案2】：

编写自己的函数来执行此操作实际上并不难 - 这是我的，欢迎您使用。

如果您在大量点上执行此操作并且速度很重要，我猜这个函数将在很长一段时间内击败基于 for-loop 的速度解决方案 - numpy 被设计为在执行操作时高效在整个矩阵上。

import numpy
a = numpy.array([1,1,1])
b = numpy.random.rand(20,3)

def euclidean_distances(ref_point, co_ords_array):
    diffs = co_ords_array - ref_point
    sqrd_diffs = numpy.square(diffs)
    sum_sqrd_diffs = numpy.sum(sqrd_diffs, axis = 1)
    euc_dists = numpy.sqrt(sum_sqrd_diffs)
    return euc_dists

【讨论】：

我完全支持你的解决方案，因此给了你一个赞成票，因为我也更喜欢自己编写这样的计算。如我的solution to this question 所示，写出计算效率更高（并且可以浓缩为单行）。我将其添加到我的解决方案中。我之所以使用scipy.spatial.distance 模块，是因为它的一大优势：只需更改函数调用，即可轻松灵活地使用各种空间距离计算。
我非常喜欢这个，虽然我同意@Scotty1-，但是能够使用距离库是很有价值的。我尽量不重写代码，但有时这是必要的，甚至是可取的。根据正在完成的工作，它可能会非常令人满意。谢谢。

【解决方案3】：

这段代码将获得欧几里得范数，该范数在许多情况下都应该有效，而且相当快，而且只有一行。根据需要，其他方法更有效或更灵活，我更喜欢根据正在完成的工作发布的其他一些解决方案。

import numpy
a = numpy.array([1,1,1])
b = numpy.random.rand(20,3)

distances = numpy.linalg.norm(a - b, axis = 1)

【讨论】：

【解决方案4】：

注意a定义中额外的[]集合

import numpy, scipy.spatial.distance
a = numpy.array([[1,1,1]])
b = numpy.random.rand(20,3)

distances = scipy.spatial.distance.cdist(b, a, metric='euclidean')

【讨论】：

聪明，你能解释一下为什么这样有效吗？谢谢。