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><strong>Approach #1</strong></p>

<p>You can simulate that iterator dependency criteria for a vectorized solution using a <code>triangular matrix</code>. This is based on <a href="http://stackoverflow.com/a/36045511/3293881"><code>this post</code></a> that dealt with multiplication involving <code>iterator dependency</code>. For performing the elementwise equality of each element in <code>vect</code> against its all elements, we can use <code>NumPy broadcasting</code>. Finally, we can use <code>np.count_nonzero</code> to get the count, as it's supposed to be very efficient in summing purposes on boolean arrays.</p>

<p>So, we would have a solution like so -</p>

<pre><code>mask = np.triu(vect[:,None] == vect,1)
counter = np.count_nonzero(mask)
dupl = np.where(mask)[1]
</code></pre>

<p>If you only care about the count <code>counter</code>, we could have two more approaches as listed next.</p>

<p><strong>Approach #2</strong></p>

<p>We can avoid the use of the triangular matrix and simply get the entire count and just subtract the contribution from diagonal elements and consider just one of either lower of upper triangular regions by just halving the remaining count as the contributions from either ones would be identical.</p>

<p>So, we would have a modified solution like so -</p>

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

<p><strong>Approach #3</strong></p>

<p>Here's an entirely different approach that uses the fact the count of each unique element plays a cumsumed contribution to the final total. </p>

<p>So, with that idea in mind, we would have a third approach like so -</p>

<pre><code>count = np.bincount(vect) # OR np.unique(vect,return_counts=True)[1]
idx = count[count&gt;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())
</code></pre>

python arrays loops numpy optimization