Vectorizing Iterative Additions to NumPy Arrays

Question

For each element in a randomized array of 2D indexes (with potential duplicates), I want "+ = 1" to the corresponding grid in a 2D zero matrix. However, I do not know how to optimize the calculations. Using the standard for the loop, as shown here,

def interadd(): U = 100 input = np.random.random(size=(5000,2)) * U idx = np.floor(input).astype(np.int) grids = np.zeros((U,U)) for i in range(len(input)): grids[idx[i,0],idx[i,1]] += 1 return grids 

The lead time can be quite significant:

 >> timeit(interadd, number=5000) 43.69953393936157 

Is there a way to vectorize this iterative process?

Answer 1

You can speed it up a bit using np.add.at , which handles the case of duplicate indexes correctly:

 def interadd(U, idx): grids = np.zeros((U,U)) for i in range(len(idx)): grids[idx[i,0],idx[i,1]] += 1 return grids def interadd2(U, idx): grids = np.zeros((U,U)) np.add.at(grids, idx.T.tolist(), 1) return grids def interadd3(U, idx): # YXD suggestion grids = np.zeros((U,U)) np.add.at(grids, (idx[:,0], idx[:,1]), 1) return grids 

which gives

 >>> U = 100 >>> idx = np.floor(np.random.random(size=(5000,2))*U).astype(np.int) >>> (interadd(U, idx) == interadd2(U, idx)).all() True >>> %timeit interadd(U, idx) 100 loops, best of 3: 8.48 ms per loop >>> %timeit interadd2(U, idx) 100 loops, best of 3: 2.62 ms per loop 

---

And the suggestion of YXD:

 >>> (interadd(U, idx) == interadd3(U, idx)).all() True >>> %timeit interadd3(U, idx) 1000 loops, best of 3: 1.09 ms per loop 

Answer 2

Divakar's answer will make me try the following, which looks like the fastest way:

 lin_idx = idx[:,0]*U + idx[:,1] grids = np.bincount(lin_idx, minlength=U**2).reshape(U, U) 

Timings:

 In [184]: U = 100 ...: input = np.random.random(size=(5000,2)) * U ...: idx = np.floor(input).astype(np.int) In [185]: %timeit interadd3(U, idx) # By DSM / XYD 1000 loops, best of 3: 1.68 ms per loop In [186]: %timeit unique_counts(U, idx) # By Divakar 1000 loops, best of 3: 676 µs per loop In [187]: %%timeit ...: lin_idx = idx[:,0]*U + idx[:,1] ...: grids = np.bincount(lin_idx, minlength=U*U).reshape(U, U) ...: 10000 loops, best of 3: 97.5 µs per loop 

Answer 3

You can convert the R,C indices from idx to linear indices, and then find the unique ones along with their counts and finally save them in grids output as the final output. Here's the implementation to achieve the same -

 # Calculate linear indices corressponding to idx lin_idx = idx[:,0]*U + idx[:,1] # Get unique linear indices and their counts unq_lin_idx,idx_counts = np.unique(lin_idx,return_counts=True) # Setup output array and store index counts in raveled/flattened version grids = np.zeros((U,U)) grids.ravel()[unq_lin_idx] = idx_counts 

Runtime Tests -

The following are run-time tests covering all approaches (including @DSM suitable ) and using the same definitions that are specified in this solution -

 In [63]: U = 100 ...: idx = np.floor(np.random.random(size=(5000,2))*U).astype(np.int) ...: In [64]: %timeit interadd(U, idx) 100 loops, best of 3: 7.57 ms per loop In [65]: %timeit interadd2(U, idx) 100 loops, best of 3: 2.59 ms per loop In [66]: %timeit interadd3(U, idx) 1000 loops, best of 3: 1.24 ms per loop In [67]: def unique_counts(U, idx): ...: lin_idx = idx[:,0]*U + idx[:,1] ...: unq_lin_idx,idx_counts = np.unique(lin_idx,return_counts=True) ...: grids = np.zeros((U,U)) ...: grids.ravel()[unq_lin_idx] = idx_counts ...: return grids ...: In [68]: %timeit unique_counts(U, idx) 1000 loops, best of 3: 595 µs per loop 

Time series show that the proposed np.unique based np.unique more than 50% faster than the second fastest approach.