numpy配列のエントリを互いに効率的に比較する方法は？

Question

Iは、各エントリが3Dベクトルである長さtot_vecのnumpyのアレイembed_vec有する：この配列の各要素について

[[ 0.52483319 0.78015841 0.71117216] 
[ 0.53041481 0.79462171 0.67234534] 
[ 0.53645428 0.80896727 0.63119403] 
..., 
[ 0.72283509 0.40070804 0.15220522] 
[ 0.71277758 0.38498613 0.16141834] 
[ 0.70221445 0.36918032 0.17370776]] 

、私が「近くにある他のエントリの数を知りたいと"そのエントリに。閉じるとは、2つのベクトルの間の距離が指定された値Rより小さいことを意味します。このためには、この配列内のすべての可能なペアを互いに比較し、配列内の各ベクトルの近接ベクトルの数を調べる必要があります。だから私はこれをやっている：私は2つのネストされたPythonのループを持っていると、より大きなアレイサイズのために、これは永遠にかかるため

p = np.zeros(tot_vec) # This contains the number of close vectors 
for i in range(tot_vec-1): 
    for j in range(i+1, tot_vec): 
     if np.linalg.norm(embed_vec[i]-embed_vec[j]) < R: 
      p[i] += 1 

しかし、これは非常に非効率的です。これがC++やFortranであれば大きな問題にはなりませんでした。私の質問は、numpyを効率的に使用するためにベクトル化手法を使って効率的に同じことを実現できますか？副次的なこととして、私はパンダを使った解決策も気にしない。

Answer 1

アプローチ＃1：ベクトル化アプローチ -

def vectorized_app(embed_vec, R): 
    tot_vec = embed_vec.shape[0]   
    r,c = np.triu_indices(tot_vec,1) 
    subs = embed_vec[r] - embed_vec[c] 
    dists = np.einsum('ij,ij->i',subs,subs) 
    return np.bincount(r,dists<R**2,minlength=tot_vec) 

アプローチ＃2：以下ループ複雑で（非常に大きなアレイの場合） -

def loopy_less_app(embed_vec, R): 
    tot_vec = embed_vec.shape[0] 
    Rsq = R**2 
    out = np.zeros(tot_vec,dtype=int) 
    for i in range(tot_vec): 
     subs = embed_vec[i] - embed_vec[i+1:tot_vec] 
     dists = np.einsum('ij,ij->i',subs,subs) 
     out[i] = np.count_nonzero(dists < Rsq) 
    return out 

---

ベンチマーク

Ori ginalアプローチ -

def loopy_app(embed_vec, R): 
    tot_vec = embed_vec.shape[0] 
    p = np.zeros(tot_vec) # This contains the number of close vectors 
    for i in range(tot_vec-1): 
     for j in range(i+1, tot_vec): 
      if np.linalg.norm(embed_vec[i]-embed_vec[j]) < R: 
       p[i] += 1 
    return p     

タイミング - そこ

In [76]: # Sample random array 
    ...: embed_vec = np.random.rand(3000,3) 
    ...: R = 0.5 
    ...: 

In [77]: %timeit loopy_app(embed_vec, R) 
1 loops, best of 3: 50.5 s per loop 

In [78]: %timeit loopy_less_app(embed_vec, R) 
10 loops, best of 3: 143 ms per loop 

350x+スピードアップ！データセットがどのように大きなに応じて、あなたがしたいことがあり、今

disp_vecs=tot_vec[:,None,:]-tot_vec[None,:,:]

：

In [81]: # Sample random array 
    ...: embed_vec = np.random.rand(20000,3) 
    ...: R = 0.5 
    ...: 

In [82]: %timeit loopy_less_app(embed_vec, R) 
1 loops, best of 3: 4.47 s per loop 

Answer 2

第1の差分を放送 - 提案loopy_less_appはとはるかに大きな配列で行く

すべての数学を除いて、拳のパス。距離がr未満の場合、すべてのコンポーネントが

disps=np.linlg.norm(disp_vec[first_mask],axis=-1) 
second_mask=disps<r 

今すぐ

disps=disps[second_mask] 
first_mask[first_mask]=second_mask 

dispsは今再割り当て未満r

first_mask=np.max(disp_vec, axis=-1)<r

は次に実際の計算を行うべきです良い値であり、first_maskはwherのブール値マスクです彼らは行く。そこから処理することができます。

Answer 3

私はこの質問に興味を持ち、scipyのcKDTreeを使って効率的に解決しようとしました。しかし、距離が< = Rのすべてのペアのリストが内部的に保持されているため、このアプローチではメモリが不足する可能性があります。あなたのRとtot_vecが十分に小さい場合、それは動作します：ケース・メモリに

import numpy as np 
from scipy.spatial import cKDTree as KDTree 

tot_vec = 60000 
embed_vec = np.random.randn(tot_vec, 3) 
R = 0.1 

tree = KDTree(embed_vec, leafsize=100) 
p = np.zeros(tot_vec) 
for pair in tree.query_pairs(R): 
    p[pair[0]] += 1 
    p[pair[1]] += 1 

が問題である、いくつかの努力でCのパフォーマンスを犠牲にPythonでジェネレータ関数としてquery_pairsを書き換えることが可能です。