Vector operations with numpy

Question

I have three numpy arrays:

X: a 3073 x 49000 matrix
W: a 10 x 3073 matrix
y: a 49000 x 1 vector

y contains values between 0 and 9, each value represents a row in W.

I would like to add the first column of X to the row in W given by the first element in y. I.e. if the first element in y is 3, add the first column of X to the fourth row of W. And then add the second column of X to the row in W given by the second element in y and so on, until all columns of X has been aded to the row in W specified by y, which means a total of 49000 added rows.

W[y] += X.T does not work for me, because this will not add more than one vector to a row in W.

Please note: I'm only looking for vectorized solutions. I.e. no for-loops.

EDIT: To clarify I'll add an example with small matrix sizes adapted from Salvador Dali's example below.

In [1]: import numpy as np

In [2]: a, b, c = 3, 4, 5

In [3]: np.random.seed(0)

In [4]: X = np.random.randint(10, size=(b,c))

In [5]: W = np.random.randint(10, size=(a,b))

In [6]: y = np.random.randint(a, size=(c,1))

In [7]: X
Out[7]: 
array([[5, 0, 3, 3, 7],
       [9, 3, 5, 2, 4],
       [7, 6, 8, 8, 1],
       [6, 7, 7, 8, 1]])

In [8]: W
Out[8]: 
array([[5, 9, 8, 9],
       [4, 3, 0, 3],
       [5, 0, 2, 3]])

In [9]: y
Out[9]: 
array([[0],
       [1],
       [1],
       [2],
       [0]])

In [10]: W[y.ravel()] + X.T
Out[10]: 
array([[10, 18, 15, 15],
       [ 4,  6,  6, 10],
       [ 7,  8,  8, 10],
       [ 8,  2, 10, 11],
       [12, 13,  9, 10]])

In [11]: W[y.ravel()] = W[y.ravel()] + X.T

In [12]: W
Out[12]: 
array([[12, 13,  9, 10],
       [ 7,  8,  8, 10],
       [ 8,  2, 10, 11]])

The problem is to get BOTH column 0 and column 4 in X added to row 0 in W, as well as both column 1 and 2 in X added to row 1 in W.

The desired outcome is thus:

W = [[17, 22, 16, 16],
     [ 7, 11, 14, 17],
     [ 8,  2, 10, 11]]

Answer 1

Vectorized approaches

Approach #1

Based on this answer, here's a vectorized solution using np.bincount -

N = y.max()+1
id = y.ravel() + np.arange(X.shape[0])[:,None]*N
W[:N] += np.bincount(id.ravel(), weights=X.ravel()).reshape(-1,N).T

Approach #2

You can make good usage of boolean indexing and np.einsum to get the job done in a concise vectorized manner -

N = y.max()+1
W[:N] += np.einsum('ijk,lk->il',(np.arange(N)[:,None,None] == y.ravel()),X)

---

Loopy approaches

Approach #3

Since you are selecting and adding up a huge number of columns from X per unique y, it might be better in terms of performance to run a loop with complexity equal to the number of such unique y's, which seems to be at max equal to the number of rows in W and that in your case is just 10. Thus, the loop has just 10 iterations, not bad! Here's the implementation to fulfill those aspirations -

for k in range(W.shape[0]):
    W[k] += X[:,(y==k).ravel()].sum(1)

Approach #4

You can bring in np.einsum to do the columnwise summations and have the final output like so -

for k in range(W.shape[0]):
    W[k] += np.einsum('ij->i',X[:,(y==k).ravel()])

Answer 2

First the straight forward loop solution as reference:

In [65]: for i,j in enumerate(y):
    W[j]+=X[:,i]
   ....:     

In [66]: W
Out[66]: 
array([[17, 22, 16, 16],
       [ 7, 11, 14, 17],
       [ 8,  2, 10, 11]])

An add.at solution:

In [67]: W=W1.copy()
In [68]: np.add.at(W,(y.ravel()),X.T)
In [69]: W
Out[69]: 
array([[17, 22, 16, 16],
       [ 7, 11, 14, 17],
       [ 8,  2, 10, 11]])

add.at does an unbuffered calculation, getting around the buffering that prevents W[y.ravel()] += X.T from working. It is still iterative, but the loop has been moved to compiled code. It isn't true vectorization because the order of application matters. The addition for one row of X.T depends on the results from the previous rows.

https://stackoverflow.com/a/20811014/901925 is the answer I gave a couple of years ago to a similar question (for 1d arrays).

But when dealing with your large arrays:

X: a 3073 x 49000 matrix
W: a 10 x 3073 matrix
y: a 49000 x 1 vector 

this can run into speed issues. Note that W[y.ravel()] is the same size as X.T (why did you pick these sizes that require transpose?). And it's a copy, not a view. So there's already a time penalty.

bincount has been suggested in previous questions, and I think it is faster. Making for loop with index arrays faster (both bincount and add.at solutions)

Iterating over the small 3073 dimension could also have speed advantages. Or better yet on the size 10 dimension as Divakar demonstrates.

---

For the small test case, a,b,c=3,4,5, the add.at solution is fastest, with Divakar's bincount and einseum next. For a larger a,b,c=10,1000,20000, add.at gets very slow, with bincount being the fastest.

Related SO answers

https://stackoverflow.com/a/28205888/901925 (notes that bincount requires complete coverage for y).

https://stackoverflow.com/a/30041823/901925 (where Divakar again shows that bincount rules!)

Answer 3

This will achieve what you want: X + W[y.ravel()].T

---

To see that this really does the work, here is a reproducible example:

import numpy as np
np.random.seed(0)
a, b, c = 3, 5, 4  # you can use your 3073, 49000, 10 later

X = np.random.rand(a, b)
W = np.random.rand(c, a)
y = np.random.randint(c, size=(b, 1))

Now your matrices are:

[[ 0.0871293   0.0202184   0.83261985]
 [ 0.77815675  0.87001215  0.97861834]
 [ 0.79915856  0.46147936  0.78052918]
 [ 0.11827443  0.63992102  0.14335329]]

[[3]
 [0]
 [3]
 [2]
 [0]]

[[ 0.5488135   0.71518937  0.60276338  0.54488318  0.4236548 ]
 [ 0.64589411  0.43758721  0.891773    0.96366276  0.38344152]
 [ 0.79172504  0.52889492  0.56804456  0.92559664  0.07103606]]

And W[y.ravel()] gives you " W given by the first element in y". By transposing it, you will get a matrix ready to be added to X:

[[ 0.11827443  0.0871293   0.11827443  0.79915856  0.0871293 ]
 [ 0.63992102  0.0202184   0.63992102  0.46147936  0.0202184 ]
 [ 0.14335329  0.83261985  0.14335329  0.78052918  0.83261985]]

Answer 4

While I can't say that this is very pythonic, it is a solution (I think):

for column in range(x.shape[1]):
    w[y[column]] = x[:,column].T