Geek Answers Handbook

Replace for numpy cast with scipy.sparse.csc_matrix

I have the following expression in my code:

a = (b / x[:, np.newaxis]).sum(axis=1)

where b is an ndarray of form (M, N) and x is an ndarray of form (M,) . Now b is actually sparse, so for memory efficiency I would like to replace with scipy.sparse.csc_matrix or csr_matrix . However, translation is not implemented in this way (even if division or multiplication guarantees the preservation of sparseness) (records x are nonzero) and raises a NotImplementedError . Is there a sparse function that I don't know that will do what I want? ( dot() will be summed on the wrong axis.)

+6
source share
2 answers

If b is in CSC format, then b.data has nonzero entries b , and b.indices has a row index for each nonzero entry, so you can do the splitting as:

 b.data /= np.take(x, b.indices)

This is a different hacker solution than Warren, but it will probably be faster in most settings:

 b = sps.rand(1000, 1000, density=0.01, format='csc') x = np.random.rand(1000) def row_divide_col_reduce(b, x): data = b.data.copy() / np.take(x, b.indices) ret = sps.csc_matrix((data, b.indices.copy(), b.indptr.copy()), shape=b.shape) return ret.sum(axis=1) def row_divide_col_reduce_bis(b, x): d = sps.spdiags(1.0/x, 0, len(x), len(x)) return (d * b).sum(axis=1) In [2]: %timeit row_divide_col_reduce(b, x) 1000 loops, best of 3: 210 us per loop In [3]: %timeit row_divide_col_reduce_bis(b, x) 1000 loops, best of 3: 697 us per loop In [4]: np.allclose(row_divide_col_reduce(b, x), ...: row_divide_col_reduce_bis(b, x)) Out[4]: True

You can reduce the time by almost half in the above example if you are doing the division in place, i.e.:

 def row_divide_col_reduce(b, x): b.data /= np.take(x, b.indices) return b.sum(axis=1) In [2]: %timeit row_divide_col_reduce(b, x) 10000 loops, best of 3: 131 us per loop
+5
source

To implement a = (b / x[:, np.newaxis]).sum(axis=1) , you can use a = b.sum(axis=1).A1 / x . The A1 attribute returns 1D ndarray, so the result is 1D ndarray, not matrix . This squeezed expression works because you both scale on x and sum along axis 1. For example:

 In [190]: b Out[190]: <3x3 sparse matrix of type '<type 'numpy.float64'>' with 5 stored elements in Compressed Sparse Row format> In [191]: bA Out[191]: array([[ 1., 0., 2.], [ 0., 3., 0.], [ 4., 0., 5.]]) In [192]: x Out[192]: array([ 2., 3., 4.]) In [193]: b.sum(axis=1).A1 / x Out[193]: array([ 1.5 , 1. , 2.25])

More generally, if you want to scale the rows of a sparse matrix with the vector x , you can multiply b left with a sparse matrix containing 1.0/x on the diagonal. The scipy.sparse.spdiags function can be used to create such a matrix. For instance:

 In [71]: from scipy.sparse import csc_matrix, spdiags In [72]: b = csc_matrix([[1,0,2],[0,3,0],[4,0,5]], dtype=np.float64) In [73]: bA Out[73]: array([[ 1., 0., 2.], [ 0., 3., 0.], [ 4., 0., 5.]]) In [74]: x = array([2., 3., 4.]) In [75]: d = spdiags(1.0/x, 0, len(x), len(x)) In [76]: dA Out[76]: array([[ 0.5 , 0. , 0. ], [ 0. , 0.33333333, 0. ], [ 0. , 0. , 0.25 ]]) In [77]: p = d * b In [78]: pA Out[78]: array([[ 0.5 , 0. , 1. ], [ 0. , 1. , 0. ], [ 1. , 0. , 1.25]]) In [79]: a = p.sum(axis=1) In [80]: a Out[80]: matrix([[ 1.5 ], [ 1. ], [ 2.25]])
+4
source

More articles:


All Articles