2013-08-30 13 views
5

私は3次元座標の3xNの配列を持っており、すべてのエントリの距離行列を効率的に計算したいと思います。 ネストされたループではなく、効率的なループ戦略が適用できますか?Numpy効率的なものすべてに対して

現在の擬似コードの実装:

for i,coord in enumerate(coords): 
    for j,coords2 in enumerate(coords): 
     if i != j: 
      dist[i,j] = numpy.norm(coord - coord2) 
+3

['scipy.spatial.distance.pdist'](http://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.distance.pdist.html#scipy.spatial.distance。 pdist)。 – BrenBarn

答えて

5

正確にあなたの結果を再現するには、次の

>>> import scipy.spatial as sp 
>>> import numpy as np 
>>> a=np.random.rand(5,3) #Note this is the transpose of your array. 
>>> a 
array([[ 0.83921304, 0.72659404, 0.50434178], #0 
     [ 0.99883826, 0.91739731, 0.9435401 ], #1 
     [ 0.94327962, 0.57665875, 0.85853404], #2 
     [ 0.30053567, 0.44458829, 0.35677649], #3 
     [ 0.01345765, 0.49247883, 0.11496977]]) #4 
>>> sp.distance.cdist(a,a) 
array([[ 0.  , 0.50475862, 0.39845025, 0.62568048, 0.94249268], 
     [ 0.50475862, 0.  , 0.35554966, 1.02735895, 1.35575051], 
     [ 0.39845025, 0.35554966, 0.  , 0.82602847, 1.1935422 ], 
     [ 0.62568048, 1.02735895, 0.82602847, 0.  , 0.3783884 ], 
     [ 0.94249268, 1.35575051, 1.1935422 , 0.3783884 , 0.  ]]) 

を計算を複製することなく、より効率的にそれを行うだけのユニークなペアを計算するには:

>>> sp.distance.pdist(a) 
array([ 0.50475862, 0.39845025, 0.62568048, 0.94249268, 0.35554966, 
     1.02735895, 1.35575051, 0.82602847, 1.1935422 , 0.3783884 ]) 
#pairs: [(0, 1), (0, 2), (0, 3), (0, 4), (1, 2), (1, 3), (1, 4), (2, 3), 
#   (2, 4), (3, 4)] 

注意2つの配列の関係あなたが望むのであれば

arr=np.random.rand(1000,3) 

%timeit sp.distance.pdist(arr) 
100 loops, best of 3: 4.26 ms per loop 

%timeit sp.distance.cdist(arr,arr) 
100 loops, best of 3: 9.31 ms per loop 

%timeit pdist_toarray(arr) 
10 loops, best of 3: 66.2 ms per loop 

%timeit looping(arr) 
1 loops, best of 3: 16.7 s per loop 

>>> out=np.zeros((a.shape[0],a.shape[0])) 
>>> dists=sp.distance.pdist(a) 
>>> out[np.triu_indices(a.shape[0],1)]=dists 
>>> out+=out.T 

>>> out 
array([[ 0.  , 0.50475862, 0.39845025, 0.62568048, 0.94249268], 
     [ 0.50475862, 0.  , 0.35554966, 1.02735895, 1.35575051], 
     [ 0.39845025, 0.35554966, 0.  , 0.82602847, 1.1935422 ], 
     [ 0.62568048, 1.02735895, 0.82602847, 0.  , 0.3783884 ], 
     [ 0.94249268, 1.35575051, 1.1935422 , 0.3783884 , 0.  ]]) 

一部やや意外timings-

セットアップ:

def pdist_toarray(a): 
    out=np.zeros((a.shape[0],a.shape[0])) 
    dists=sp.distance.pdist(a) 

    out[np.triu_indices(a.shape[0],1)]=dists 
    return out+out.T 

def looping(a): 
    out=np.zeros((a.shape[0],a.shape[0])) 
    for i in xrange(a.shape[0]): 
     for j in xrange(a.shape[0]): 
      out[i,j]=np.linalg.norm(a[i]-a[j]) 
    return out 

タイミングをcdist配列を再生することができます正方形あなたがちょうどペアを使用する場合はを使用する場合は、cdistを使用する必要があります。ルーピングは、要素が1000の配列では約4000倍遅く、要素が10の配列では約70倍遅くなるのはcdistです。

+0

Scipyはベクトルと行列を変換する['squareform'](https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.spatial.distance.squareform.html)関数を提供しています距離配列のバージョン(「凝縮」および「冗長」形式)は、 – Praveen

関連する問題