numpyの効率的なループ

Question

この質問のバージョンは既に尋ねられていますが、私は満足のいく回答は見つかりませんでした。

問題：大きなnumpyベクトルがある場合、複製されるベクトル要素のインデックスを見つけます（そのバリエーションは許容誤差との比較が可能です）。

したがって、問題は〜O（N^2）とメモリの境界です（少なくとも現在のアルゴリズムの観点から）。なぜ私がPythonを試してみたとしても、同等のCコードよりも100倍以上遅いのではないかと思います。

import numpy as np 
N = 10000 
vect = np.arange(float(N)) 
vect[N/2] = 1 
vect[N/4] = 1 
dupl = [] 
print("init done") 
counter = 0 
for i in range(N): 
    for j in range(i+1, N): 
     if vect[i] == vect[j]: 
      dupl.append(j) 
      counter += 1 

print("counter =", counter) 
print(dupl) 
# For simplicity, this code ignores repeated indices 
# which can be trimmed later. Ref output is 
# counter = 3 
# [2500, 5000, 5000] 

私はnumpyの反復子を使用して試みたが、彼らはさらに悪化している（〜x4-5） http://docs.scipy.org/doc/numpy/reference/arrays.nditer.html

（上記コード）N = Iは、PythonでCで0.1秒、12秒を取得してい10,000使用np.nditerを使ってPythonで40秒、np.ndindexを使ってPythonで50秒。私はそれをN = 160,000にプッシュし、タイミングは予想どおりN^2になりました。

Answer 1

アプローチ＃1

あなたはtriangular matrixを使用してベクトル化液とイテレータ依存条件をシミュレートすることができます。これは、iterator dependencyを含む乗算を扱ったthis postに基づいています。 vectの各要素の要素ごとの等価性をすべての要素に対して実行するには、NumPy broadcastingを使用できます。最後に、ブール値配列での集計を非常に効率的にするため、np.count_nonzeroを使用してカウントを取得できます。あなただけの数counterを気にしている場合、次の記載されているとして、我々は2つのアプローチを持っている可能性があり、

mask = np.triu(vect[:,None] == vect,1) 
counter = np.count_nonzero(mask) 
dupl = np.where(mask)[1] 

-

そこで、我々はそうのような解決策を持っているでしょう。

私たちは、三角行列の使用を避け、単に全体のカウントを取得し、ちょうど対角要素からの寄与を差し引くと、ちょうど半分にすることにより、上三角部分のどちらか下のちょうど1を考慮することができる2

アプローチ＃いずれかからの寄与としての残りのカウントは同一である。

だから、私たちはそうのような修正ソリューションだろう -

counter = (np.count_nonzero(vect[:,None] == vect) - vect.size)//2 

アプローチ＃3

をここでは、各ユニークな要素のカウントがcumsumed貢献を果たしているという事実を使用しています全く異なるアプローチがあります最終的な合計に。

だから、このことを念頭に置いてアイデアを、私たちはそうのような第三のアプローチを持っているでしょう -

count = np.bincount(vect) # OR np.unique(vect,return_counts=True)[1] 
idx = count[count>1] 
id_arr = np.ones(idx.sum(),dtype=int) 
id_arr[0] = 0 
id_arr[idx[:-1].cumsum()] = -idx[:-1]+1 
counter = np.sum(id_arr.cumsum()) 

Answer 2

のPython自体は非常にダイナミック、遅い、言語です。 numpyのアイデアは、vectorizationを使用し、明示的なループを避けることです。この場合、np.equal.outerを使用できます。あなたは合計を見つけるために、例えば、今

a = np.equal.outer(vect, vect) 

で始めることができます。

>>> np.sum(a) 
10006 

等しい私のインデックスを見つけるには、

np.fill_diagonal(a, 0) 

>>> np.nonzero(np.any(a, axis=0))[0] 
array([ 1, 2500, 5000]) 

---

を行うことができます

タイミング

def find_vec(): 
    a = np.equal.outer(vect, vect) 
    s = np.sum(a) 
    np.fill_diagonal(a, 0) 
    return np.sum(a), np.nonzero(np.any(a, axis=0))[0] 

>>> %timeit find_vec() 
1 loops, best of 3: 214 ms per loop 

def find_loop(): 
    dupl = [] 
    counter = 0 
    for i in range(N): 
     for j in range(i+1, N): 
      if vect[i] == vect[j]: 
       dupl.append(j) 
       counter += 1 
    return dupl 

>>> % timeit find_loop() 
1 loops, best of 3: 8.51 s per loop 

Answer 3

私はPythonが同等のCコードよりも100倍以上遅いしようとした理由は何だろうか。

通常、PythonプログラムはCプログラムよりも100倍遅いためです。

Cで重要なコードパスを実装し、Python-Cバインディングを提供するか、アルゴリズムを変更することができます。配列を値からインデックスに逆転させるdictを使用してO（N）バージョンを書くことができます。

import numpy as np 
N = 10000 
vect = np.arange(float(N)) 
vect[N/2] = 1 
vect[N/4] = 1 
dupl = {} 
print("init done") 
counter = 0 
for i in range(N): 
    e = dupl.get(vect[i], None) 
    if e is None: 
     dupl[vect[i]] = [i] 
    else: 
     e.append(i) 
     counter += 1 

print("counter =", counter) 
print([(k, v) for k, v in dupl.items() if len(v) > 1]) 

編集：

あなたは腹筋とEPSに対してテストする必要がある場合（vectを[I] - vectを[J]）< EPSあなたはEPS

abs(vect[i] - vect[j]) < eps -> 
abs(vect[i] - vect[j])/eps < (eps/eps) -> 
abs(vect[i]/eps - vect[j]/eps) < 1 
int(abs(vect[i]/eps - vect[j]/eps)) = 0 

に値を正規化することができますこのよう

：

import numpy as np 
N = 10000 
vect = np.arange(float(N)) 
vect[N/2] = 1 
vect[N/4] = 1 
dupl = {} 
print("init done") 
counter = 0 
eps = 0.01 
for i in range(N): 
    k = int(vect[i]/eps) 
    e = dupl.get(k, None) 
    if e is None: 
     dupl[k] = [i] 
    else: 
     e.append(i) 
     counter += 1 

print("counter =", counter) 
print([(k, v) for k, v in dupl.items() if len(v) > 1]) 

Answer 4

numpy_indexedパッケージを使用したこのソリューションは、コンプレックスを持っていますn n Log nであり、完全にベクトル化されている。 Cのパフォーマンスとはまったく異なるわけではありません。

import numpy_indexed as npi 
dpl = np.flatnonzero(npi.multiplicity(vect) > 1) 

Answer 5

なぜこのようにしたいのですか？ NumPy配列は、不透明なデータ構造を意図しています。これは、NumPy配列がNumPyシステム内に作成され、NumPyサブシステムに送信されて結果が返されることを意味します。 NumPyは、あなたがリクエストを出して結果を出すブラックボックスでなければなりません。

上記のコードでは、NumPyのパフォーマンスが恐ろしいほど悪いと私はまったく驚いていません。

次はあなたが欲しいものを効果的にする必要があり、私は信じていますが、numpyの道に行わ：亜美Tavoryの答えの代わりとして

import numpy as np 

N = 10000 
vect = np.arange(float(N)) 
vect[N/2] = 1 
vect[N/4] = 1 

print([np.where(a == vect)[0] for a in vect][1]) 

# Delivers [1, 2500, 5000] 

Answer 6

を、あなたは重複を検出するために、コレクションのパッケージからCounterを使用することができます。私のコンピュータでは、それはもっと速いようです。以下の関数を参照してください。異なる関数もあります。

import collections 
import numpy as np 

def find_duplicates_original(x): 
    d = [] 
    for i in range(len(x)): 
     for j in range(i + 1, len(x)): 
      if x[i] == x[j]: 
       d.append(j) 
    return d 

def find_duplicates_outer(x): 
    a = np.equal.outer(x, x) 
    np.fill_diagonal(a, 0) 
    return np.flatnonzero(np.any(a, axis=0)) 

def find_duplicates_counter(x): 
    counter = collections.Counter(x) 
    values = (v for v, c in counter.items() if c > 1) 
    return {v: np.flatnonzero(x == v) for v in values} 


n = 10000 
x = np.arange(float(n)) 
x[n // 2] = 1 
x[n // 4] = 1 

>>>> find_duplicates_counter(x) 
{1.0: array([ 1, 2500, 5000], dtype=int64)} 

>>>> %timeit find_duplicates_original(x) 
1 loop, best of 3: 12 s per loop 

>>>> %timeit find_duplicates_outer(x) 
10 loops, best of 3: 84.3 ms per loop 

>>>> %timeit find_duplicates_counter(x) 
1000 loops, best of 3: 1.63 ms per loop 

Answer 7

これは、コードの18秒と比較して8ミリ秒で実行され、奇妙なライブラリは使用しません。 @ vs0のアプローチと似ていますが、私はdefaultdictが好きです。それはおよそO（N）でなければなりません。

from collections import defaultdict 
dupl = [] 
counter = 0 
indexes = defaultdict(list) 
for i, e in enumerate(vect): 
    indexes[e].append(i) 
    if len(indexes[e]) > 1: 
     dupl.append(i) 
     counter += 1 

Answer 8

回答が終わり、全く満足できるものはなかったので、私は自分の解決策を投稿します。

これは、私が最初に思ったようにネストされたループではなく、このケースではPythonを遅くする割当だと私は理解しています。ライブラリやコンパイル済みコードを使用すると、割り当ての必要性がなくなり、パフォーマンスが大幅に向上します。 0.2秒 -

from __future__ import print_function 
import numpy as np 
from numba import jit 

N = 10000 
vect = np.arange(N, dtype=np.float32) 

vect[N/2] = 1 
vect[N/4] = 1 
dupl = np.zeros(N, dtype=np.int32) 

print("init done") 
# uncomment to enable compiled function 
#@jit 
def duplicates(i, counter, dupl, vect): 
    eps = 0.01 
    ns = len(vect) 
    for j in range(i+1, ns): 
     # replace if to use approx comparison 
     #if abs(vect[i] - vect[j]) < eps: 
     if vect[i] == vect[j]: 
      dupl[counter] = j 
      counter += 1 
    return counter 

counter = 0 
for i in xrange(N): 
    counter = duplicates(i, counter, dupl, vect) 

print("counter =", counter) 
print(dupl[0:counter]) 

テスト

# no jit 
$ time python array-test-numba.py 
init done 
counter = 3 
[2500 5000 5000] 

elapsed 10.135 s 

# with jit 
$ time python array-test-numba.py 
init done 
counter = 3 
[2500 5000 5000] 

elapsed 0.480 s 

（コメントアウト@jitで）コンパイルされたバージョンの性能は、Cコードのパフォーマンス〜0.1に近いです。最終ループを削除すると、パフォーマンスがさらに向上する可能性があります。パフォーマンスの差は、コンパイルされたバージョンの違いはほとんどないが、epsを使用した近似比較を使用するとさらに強力になります。

# no jit 
$ time python array-test-numba.py 
init done 
counter = 3 
[2500 5000 5000] 

elapsed 109.218 s 

# with jit 
$ time python array-test-numba.py 
init done 
counter = 3 
[2500 5000 5000] 

elapsed 0.506 s 

これは〜200xの違いです。実際のコードでは、両方のループを関数に入れて、変数型の関数テンプレートを使用する必要がありました。少し複雑ですが、それほど多くはありませんでした。