在numpy/scipy中高效地计算N个点与参考点之间的距离

Question

我刚开始使用scipy/numpy。我有一个100000*3的数组，每一行都是一个坐标，还有一个1*3的中心点。我想计算数组中每一行到中心的距离并将它们存储在另一个数组中。最有效的方法是什么？

- D. Huang

---

可能是使用numpy计算欧几里得距离的重复问题。 - Fred Foo

@larsmans：我认为这不是重复，因为答案只涉及两点之间的距离，而不是N个点与参考点之间的距离。当然，这些回答没有指向我下面展示的高效scipy解决方案。 - JoshAdel

Answer 1

我会建议您查看 scipy.spatial.distance.cdist：

http://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.distance.cdist.html

import numpy as np
import scipy

a = np.random.normal(size=(10,3))
b = np.random.normal(size=(1,3))

dist = scipy.spatial.distance.cdist(a,b) # pick the appropriate distance metric 

dist 的默认距离度量等价于：

np.sqrt(np.sum((a-b)**2,axis=1))  

虽然对于大型数组，cdist更加高效（对于您的问题大小，对于我的机器而言，cdist快了约35倍）。

- JoshAdel

---

在这个答案中，单一的参考点在哪里？ - MikeB

b 是三维空间中的单一参考点，a 是另外10个三维空间中的点。 - Aaron Bramson

如果b有更多的点（对），则为np.sqrt(np.sum((hs[:, None] - an)**2, axis=2))。 - Ariel

Answer 2

我会使用sklearn实现欧氏距离（Euclidean distance）。优点在于使用矩阵乘法表达式，更加高效：

dist(x, y) = sqrt(np.dot(x, x) - 2 * np.dot(x, y) + np.dot(y, y)

一个简单的脚本看起来像这样：

import numpy as np

x = np.random.rand(1000, 3)
y = np.random.rand(1000, 3)

dist = np.sqrt(np.dot(x, x)) - (np.dot(x, y) + np.dot(x, y)) + np.dot(y, y)

这种方法的优点已经在sklearn文档中得到了很好的描述：http://scikit-learn.org/stable/modules/generated/sklearn.metrics.pairwise.euclidean_distances.html#sklearn.metrics.pairwise.euclidean_distances。

我正在使用这种方法来处理大型数据矩阵（10000, 10000），并进行了一些小的修改，如使用np.einsum函数。

- lemitech

---

不涉及针对单个参考点进行计算的问题。 - drewid

numpy.sqrt((X**2).sum(axis=1)[:, None] - 2 * X.dot(Y.transpose()) + ((Y**2).sum(axis=1)[None, :]) - BGabor

Answer 3

你还可以使用范数的发展（类似于显着的身份）。这可能是计算点阵距离最有效的方法。

这里是我最初在Octave中用于k-最近邻实现的代码片段，但你可以很容易地将其适应于numpy，因为它只使用矩阵乘法（相当于numpy.dot()）：

% Computing the euclidian distance between each known point (Xapp) and unknown points (Xtest)
% Note: we use the development of the norm just like a remarkable identity:
% ||x1 - x2||^2 = ||x1||^2 + ||x2||^2 - 2*<x1,x2>
[napp, d] = size(Xapp);
[ntest, d] = size(Xtest);

A = sum(Xapp.^2, 2);
A = repmat(A, 1, ntest);

B = sum(Xtest.^2, 2);
B = repmat(B', napp, 1);

C = Xapp*Xtest';

dist = A+B-2.*C;

- gaborous

Answer 4

这可能不是直接回答您问题的答案，但如果您想要所有粒子对的排列组合，我发现以下解决方案在某些情况下比pdist函数更快。

import numpy as np

L   = 100       # simulation box dimension
N   = 100       # Number of particles
dim = 2         # Dimensions

# Generate random positions of particles
r = (np.random.random(size=(N,dim))-0.5)*L

# uti is a list of two (1-D) numpy arrays  
# containing the indices of the upper triangular matrix
uti = np.triu_indices(100,k=1)        # k=1 eliminates diagonal indices

# uti[0] is i, and uti[1] is j from the previous example 
dr = r[uti[0]] - r[uti[1]]            # computes differences between particle positions
D = np.sqrt(np.sum(dr*dr, axis=1))    # computes distances; D is a 4950 x 1 np array

请参阅我的博客文章this，以深入了解此问题。

- shrokmel

Answer 5

你可能需要更详细地指定你感兴趣的距离函数的方式，但以下是一个非常简单（且高效）的欧几里得距离平方实现方法，基于内积 （显然可以直接推广到其他类型的距离度量）：

In []: P, c= randn(5, 3), randn(1, 3)
In []: dot(((P- c)** 2), ones(3))
Out[]: array([  8.80512,   4.61693,   2.6002,   3.3293,  12.41800])

其中P代表你的点，c代表中心。

- eat

---

在我的电脑上，这仍然比OP问题规模的“cdist”慢18倍。 - JoshAdel

@JoshAdel：这是很大的差异。顺便说一下，在我的普通机器上，使用numpy 1.6，当n= 1e5时，计时为cdist 3.5毫秒和dot 9.5毫秒。因此，dot只慢了约3倍。但是，对于更小的n（<2e3），'dot'将更快。谢谢。 - eat

Answer 6

#is it true, to find the biggest distance between the points in surface?

from math import sqrt

n = int(input( "enter the range : "))
x = list(map(float,input("type x coordinates: ").split()))
y = list(map(float,input("type y coordinates: ").split()))
maxdis = 0  
for i in range(n):
    for j in range(n):
        print(i, j, x[i], x[j], y[i], y[j])
        dist = sqrt((x[j]-x[i])**2+(y[j]-y[i])**2)
        if maxdis < dist:

            maxdis = dist
print(" maximum distance is : {:5g}".format(maxdis))

- Mohammad Feroz

---

请解释一下你的解决方案。 - Lithilion