efficient loop over numpy array
<p>Versions of this question have already been asked but I have not found a satisfactory answer.</p>

<p><strong>Problem</strong>: given a large numpy vector, find indices of the vector elements which are duplicated (a variation of that could be comparison with tolerance). </p>

<p>So the problem is ~O(N^2) and memory bound (at least from the current algorithm point of view). I wonder why whatever I tried Python is 100x or more slower than an equivalent C code.</p>

<pre><code>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]
</code></pre>

<p>I tried using numpy iterators but they are even worse (~ x4-5)
<a href="http://docs.scipy.org/doc/numpy/reference/arrays.nditer.html" rel="nofollow">http://docs.scipy.org/doc/numpy/reference/arrays.nditer.html</a></p>

<p>Using N=10,000 I'm getting 0.1 sec in C, 12 sec in Python (code above), 40 sec in Python using np.nditer, 50 sec in Python using np.ndindex. I pushed it to N=160,000 and the timing scales as N^2 as expected.</p>

<p>Python itself is a highly-dynamic, slow, language. The idea in numpy is to use <a href="https://www.safaribooksonline.com/library/view/python-for-data/9781449323592/ch04.html" rel="nofollow">vectorization</a>, and avoid explicit loops. In this case, you can use <a href="http://docs.scipy.org/doc/numpy/reference/generated/numpy.ufunc.outer.html" rel="nofollow"><code>np.equal.outer</code></a>. You can start with</p>

<pre><code>a = np.equal.outer(vect, vect)
</code></pre>

<p>Now, for example, to find the sum:</p>

<pre><code> &gt;&gt;&gt; np.sum(a)
 10006
</code></pre>

<p>To find the indices of <em>i</em> that are equal, you can do</p>

<pre><code>np.fill_diagonal(a, 0)

&gt;&gt;&gt; np.nonzero(np.any(a, axis=0))[0]
array([   1, 2500, 5000])
</code></pre>

<hr>

<p><strong>Timing</strong></p>

<pre><code>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]

&gt;&gt;&gt; %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

&gt;&gt;&gt; % timeit find_loop()
1 loops, best of 3: 8.51 s per loop
</code></pre>

python arrays loops numpy optimization