Scipy Sparse Cumsum

Question

Suppose I have a scipy.sparse.csr_matrix representing the values below

[[0 0 1 2 0 3 0 4]
 [1 0 0 2 0 3 4 0]]

I want to calculate the cumulative sum of non-zero values in-place, which would change the array to:

[[0 0 1 3 0 6 0 10]
 [1 0 0 3 0 6 10 0]]

The actual values are not 1, 2, 3, ...

The number of non-zero values in each row are unlikely to be the same.

How to do this fast?

Current program:

import scipy.sparse
import numpy as np

# sparse data
a = scipy.sparse.csr_matrix(
    [[0,0,1,2,0,3,0,4],
     [1,0,0,2,0,3,4,0]], 
    dtype=int)

# method
indptr = a.indptr
data = a.data
for i in range(a.shape[0]):
    st = indptr[i]
    en = indptr[i + 1]
    np.cumsum(data[st:en], out=data[st:en])

# print result
print(a.todense())

Result:

[[ 0  0  1  3  0  6  0 10]
 [ 1  0  0  3  0  6 10  0]]

Luthern

For working code, you should post to codereview.stackexchange.com – Hodosh

There are a lot more numpy/scipy eyes on SO than on CR. Speed questions on working code are answered all the time on SO, especially with the code packages are somewhat specialized. – Subminiature

@r xu, what you show looks good. Applying cumsum row by row is really only way to go. And your use of out is clever. There is a as strided based indptr iterator that might improve speed a bit. – Subminiature

Answer 1

How about doing this instead

a = np.array([[0,0,1,2,0,3,0,4],
              [1,0,0,2,0,3,4,0]], dtype=int)

b = a.copy()
b[b > 0] = 1
z = np.cumsum(a,axis=1)
print(z*b)

Yields

array([[ 0,  0,  1,  3,  0,  6,  0, 10],
   [ 1,  0,  0,  3,  0,  6, 10,  0]])

Doing sparse

def sparse(a):
    a = scipy.sparse.csr_matrix(a)

    indptr = a.indptr
    data = a.data
    for i in range(a.shape[0]):
        st = indptr[i]
        en = indptr[i + 1]
        np.cumsum(data[st:en], out=data[st:en])


In[1]: %timeit sparse(a)
10000 loops, best of 3: 167 µs per loop

Using multiplication

def mult(a):
    b = a.copy()
    b[b > 0] = 1
    z = np.cumsum(a, axis=1)
    z * b

In[2]: %timeit mult(a)
100000 loops, best of 3: 5.93 µs per loop

Politico

Sparse is the same code you provided and the multiplication is the code I provided – Politico

That means... a method that works on a sparse matrix, instead of a dense matrix. Your mult function currently works on a dense matrix. – Luthern