Scipy Sparse Matrix special substraction

Question

I'm doing a project and I'm doing a lot of matrix computation in it.

I'm looking for a smart way to speed up my code. In my project, I'm dealing with a sparse matrix of size 100Mx1M with around 10M non-zeros values. The example below is just to see my point.

Let's say I have:

• A vector v of size (2)
• A vector c of size (3)

• A sparse matrix X of size (2,3)

  v = np.asarray([10, 20])
  c = np.asarray([ 2,  3,  4])
  data = np.array([1, 1, 1, 1])
  row  = np.array([0, 0, 1, 1])
  col  = np.array([1, 2, 0, 2])
  X = coo_matrix((data,(row,col)), shape=(2,3))
  X.todense()
  # matrix([[0, 1, 1],
  #         [1, 0, 1]])
  

Currently I'm doing:

result = np.zeros_like(v)
d = scipy.sparse.lil_matrix((v.shape[0], v.shape[0]))
d.setdiag(v)
tmp = d * X

print tmp.todense()
#matrix([[  0.,  10.,  10.],
#        [ 20.,   0.,  20.]])
# At this point tmp is csr sparse matrix

for i in range(tmp.shape[0]):
     x_i = tmp.getrow(i)
     result += x_i.data * ( c[x_i.indices] - x_i.data)
     # I only want to do the subtraction on non-zero elements    

print result
# array([-430, -380])

And my problem is the for loop and especially the subtraction. I would like to find a way to vectorize this operation by subtracting only on the non-zero elements.

Something to get directly the sparse matrix on the subtraction:

matrix([[  0.,  -7.,  -6.],
       [ -18.,   0.,  -16.]])

Is there a way to do this smartly ?

Answer 1

You don't need to loop over the rows to do what you are already doing. And you can use a similar trick to perform the multiplication of the rows by the first vector:

import scipy.sparse as sps

# number of nonzero entries per row of X
nnz_per_row = np.diff(X.indptr)
# multiply every row by the corresponding entry of v
# You could do this in-place as:
# X.data *= np.repeat(v, nnz_per_row)
Y = sps.csr_matrix((X.data * np.repeat(v, nnz_per_row), X.indices, X.indptr),
                   shape=X.shape)

# subtract from the non-zero entries the corresponding column value in c...
Y.data -= np.take(c, Y.indices)
# ...and multiply by -1 to get the value you are after
Y.data *= -1

To see that it works, set up some dummy data

rows, cols = 3, 5
v = np.random.rand(rows)
c = np.random.rand(cols)
X = sps.rand(rows, cols, density=0.5, format='csr')

and after run the code above:

>>> x = X.toarray()
>>> mask = x == 0
>>> x *= v[:, np.newaxis]
>>> x = c - x
>>> x[mask] = 0
>>> x
array([[ 0.79935123,  0.        ,  0.        , -0.0097763 ,  0.59901243],
       [ 0.7522559 ,  0.        ,  0.67510109,  0.        ,  0.36240006],
       [ 0.        ,  0.        ,  0.72370725,  0.        ,  0.        ]])
>>> Y.toarray()
array([[ 0.79935123,  0.        ,  0.        , -0.0097763 ,  0.59901243],
       [ 0.7522559 ,  0.        ,  0.67510109,  0.        ,  0.36240006],
       [ 0.        ,  0.        ,  0.72370725,  0.        ,  0.        ]])

The way you are accumulating your result requires that there are the same number of non-zero entries in every row, which seems a pretty weird thing to do. Are you sure that is what you are after? If that's really what you want you could get that value with something like:

result = np.sum(Y.data.reshape(Y.shape[0], -1), axis=0)

but I have trouble believing that is really what you are after...