【问题标题】:Efficient distance calculation between N points and a reference in numpy/scipyN点与numpy/scipy中的参考之间的有效距离计算
【发布时间】:2011-09-19 18:44:32
【问题描述】:

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

【问题讨论】:

  • @larsmans:我不认为这是重复的,因为答案只涉及两点之间的距离,而不是 N 点和参考点之间的距离。当然,响应并没有将 OP 指向我在下面展示的有效 scipy 解决方案。

标签: python arrays numpy scipy scipy-spatial


【解决方案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 倍)。

【讨论】:

  • 在这个答案中,单一参考点在哪里?
  • b是三个维度上的单个参考点,a是三个维度上的其他10个点。
【解决方案2】:

我会使用欧几里得距离的 sklearn 实现。优点是通过使用矩阵乘法来使用更有效的表达式:

dist(x, y) = sqrt(dot(x, x) - 2 * dot(x, y) + 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)) - (dot(x, y) + dot(x, y)) + 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 函数。

【讨论】:

  • 没有解决针对单个参考点进行计算的问题
  • numpy.sqrt((X**2).sum(axis=1)[:, None] - 2 * X.dot(Y.transpose()) + ((Y**2).sum(axis=1)[None, :])
【解决方案3】:

您还可以使用发展规范(类似于显着身份)。这可能是计算点矩阵距离的最有效方法。

这是一个代码 sn-p,我最初在 Octave 中用于 k-Nearest-Neighbors 实现,但您可以轻松地将其调整为 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;

【讨论】:

    【解决方案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,在我的博文中更深入地了解这个问题。

    【讨论】:

      【解决方案5】:

      您可能需要以更详细的方式指定您感兴趣的距离函数,但这里是基于inner productSquared Euclidean Distance 的一个非常简单(高效)的实现(显然可以概括,直接的方式,到其他类型的距离度量):

      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 是中心。

      【讨论】:

      • 在我的机器上,对于 OP 的问题大小,这仍然比 cdist 慢 18 倍。
      • @JoshAdel:差别很大。 FWIW,在我的普通机器中使用numpy 1.6:对于n= 1e5,时间为cdist 3.5 ms 和dot 9.5 ms。所以 dot 只慢了大约 3 倍。然而,使用更小的n (
      【解决方案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))
      

      【讨论】:

      • 请解释您的解决方案
      猜你喜欢
      • 1970-01-01
      • 2022-01-08
      • 2023-03-23
      • 2015-11-01
      • 2015-10-29
      • 1970-01-01
      • 2014-11-29
      • 2017-04-21
      • 1970-01-01
      相关资源
      最近更新 更多