numpy: efficiently add rows of a matrix
<p>I have a matrix.  </p>

<pre><code>mat = array([
   [ 0,  1,  2,  3],
   [ 4,  5,  6,  7],
   [ 8,  9, 10, 11]
   ])
</code></pre>

<p>I'd like to get the sum of the rows at certain indices: eg.</p>

<pre><code>ixs = np.array([0,2,0,0,0,1,1])
</code></pre>

<p>I know I can compute the answer as:</p>

<pre><code>mat[ixs].sum(axis=0)
&gt; array([16, 23, 30, 37])
</code></pre>

<p>The problem is ixs may be very long, and I don't want to use all the memory to create the intermediate product mat[ixs], only to reduce it again with the sum.</p>

<p>I also know I could simply count up the indices and use multiplication instead. </p>

<pre><code>np.bincount(ixs, minlength=mat.shape[0).dot(mat)
&gt; array([16, 23, 30, 37])
</code></pre>

<p>But that will be expensive if my ixs are sparse. </p>

<p>I know about scipy's sparse matrices, and I suppose I could use them, but I'd prefer a pure numpy solution as sparse matrices are limited in various ways (such as only being 2-d)</p>

<p>So, is there a pure numpy way to merge the indexing and sum-reduction in this case?</p>

<h1>Conclusions:</h1>

<p>Thanks you Divakar and hpaulj for your very thorough responses.  By "sparse" I meant that most of the values in <code>range(w.shape[0])</code> are not in ixs.  Using that new definition (and with more realisitic data size, I re-ran Divakar tests, with some new funcitona dded :</p>

<pre><code>rng = np.random.RandomState(1234)
mat = rng.randn(1000, 500)
ixs = rng.choice(rng.randint(mat.shape[0], size=mat.shape[0]/10), size=1000)

# Divakar's solutions
In[42]: %timeit org_indexing_app(mat, ixs)
1000 loops, best of 3: 1.82 ms per loop
In[43]: %timeit org_bincount_app(mat, ixs)
The slowest run took 4.07 times longer than the fastest. This could mean that an intermediate result is being cached.
10000 loops, best of 3: 177 Âµs per loop
In[44]: %timeit indexing_modified_app(mat, ixs)
1000 loops, best of 3: 1.81 ms per loop
In[45]: %timeit bincount_modified_app(mat, ixs)
1000 loops, best of 3: 258 Âµs per loop
In[46]: %timeit simply_indexing_app(mat, ixs)
1000 loops, best of 3: 1.86 ms per loop
In[47]: %timeit take_app(mat, ixs)
1000 loops, best of 3: 1.82 ms per loop
In[48]: %timeit unq_mask_einsum_app(mat, ixs)
10 loops, best of 3: 58.2 ms per loop 
# hpaulj's solutions
In[53]: %timeit hpauljs_sparse_solution(mat, ixs)
The slowest run took 9.34 times longer than the fastest. This could mean that an intermediate result is being cached.
1000 loops, best of 3: 524 Âµs per loop
%timeit hpauljs_second_sparse_solution(mat, ixs)
100 loops, best of 3: 9.91 ms per loop
# Sparse version of original bincount solution (see below):
In[60]: %timeit sparse_bincount(mat, ixs)
10000 loops, best of 3: 71.7 Âµs per loop
</code></pre>

<p><strong>The winner</strong> in this case is the sparse version of the bincount solution.</p>

<pre><code>def sparse_bincount(mat, ixs):
    x = np.bincount(ixs)
    nonzeros, = np.nonzero(x)
    x[nonzeros].dot(mat[nonzeros])
</code></pre>

<p>An alternative to <code>bincount</code> is <code>add.at</code>:</p>

<pre><code>In [193]: mat
Out[193]: 
array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11]])
In [194]: ixs
Out[194]: array([0, 2, 0, 0, 0, 1, 1])

In [195]: J = np.zeros(mat.shape[0],int)
In [196]: np.add.at(J, ixs, 1)
In [197]: J
Out[197]: array([4, 2, 1])

In [198]: np.dot(J, mat)
Out[198]: array([16, 23, 30, 37])
</code></pre>

<p>By the sparsity, you mean, I assume, that <code>ixs</code> might not include all the rows, for example, <code>ixs</code> without the 0s:</p>

<pre><code>In [199]: ixs = np.array([2,1,1])
In [200]: J=np.zeros(mat.shape[0],int)
In [201]: np.add.at(J, ixs, 1)
In [202]: J
Out[202]: array([0, 2, 1])
In [203]: np.dot(J, mat)
Out[203]: array([16, 19, 22, 25])
</code></pre>

<p><code>J</code> still has the <code>mat.shape[0]</code> shape.  But the <code>add.at</code> should scale as the length of <code>ixs</code>.  </p>

<p>A sparse solution would look something like:</p>

<p>Make a sparse matrix from <code>ixs</code> that looks like:</p>

<pre><code>In [204]: I
Out[204]: 
array([[1, 0, 1, 1, 1, 0, 0],
       [0, 0, 0, 0, 0, 1, 1],
       [0, 1, 0, 0, 0, 0, 0]])
</code></pre>

<p>sum the rows; sparse does this with matrix multiplication like:</p>

<pre><code>In [205]: np.dot(I, np.ones((7,),int))
Out[205]: array([4, 2, 1])
</code></pre>

<p>then do our dot:</p>

<pre><code>In [206]: np.dot(np.dot(I, np.ones((7,),int)), mat)
Out[206]: array([16, 23, 30, 37])
</code></pre>

<p>Or in sparse code:</p>

<pre><code>In [225]: J = sparse.coo_matrix((np.ones_like(ixs,int),(np.arange(ixs.shape[0]), ixs)))
In [226]: J.A
Out[226]: 
array([[1, 0, 0],
       [0, 0, 1],
       [1, 0, 0],
       [1, 0, 0],
       [1, 0, 0],
       [0, 1, 0],
       [0, 1, 0]])
In [227]: J.sum(axis=0)*mat
Out[227]: matrix([[16, 23, 30, 37]])
</code></pre>

<p><code>sparse</code>, when converting from <code>coo</code> to <code>csr</code> sums duplicates.  I can take advantage that with</p>

<pre><code>In [229]: J = sparse.coo_matrix((np.ones_like(ixs,int), (np.zeros_like(ixs,int), ixs)))
In [230]: J
Out[230]: 
&lt;1x3 sparse matrix of type '&lt;class 'numpy.int32'&gt;'
    with 7 stored elements in COOrdinate format&gt;
In [231]: J.A
Out[231]: array([[4, 2, 1]])
In [232]: J*mat
Out[232]: array([[16, 23, 30, 37]], dtype=int32)
</code></pre>

python numpy indexing