I want to generate a large sparse matrix and sum it but I encounter MemoryError a lot. So I tried the operation via scipy.sparse.csc_matrix.sum instead but found that the type of data changed back to a numpy matrix after taking the sum.
window = 10 np.random.seed = 0 mat = sparse.csc_matrix(np.random.rand(100, 120)>0.5, dtype='d') print type(mat) >>> <class 'scipy.sparse.csc.csc_matrix'> mat_head = mat[:,0:window].sum(axis=1) print type(mat_head) >>> <class 'numpy.matrixlib.defmatrix.matrix'> So I generated mat as zeros matrix just to test the result when mat_head is all zeros.
mat = sparse.csc_matrix((100,120)) print type(mat) >>> <class 'scipy.sparse.csc.csc_matrix'> mat_head = mat.sum(axis=1) print type(mat_head) >>> <class 'numpy.matrixlib.defmatrix.matrix'> print np.count_nonzero(mat_head) >>> 0 Why does this happen? So sum via scipy.sparse is not benefited for preserving memory than numpy as they change the data type back anyway?
3 Answers
As far as it is possible to give a hard reason for what is essentially a design choice I'd make the following argument:
The csr and csc formats are designed for sparse but not extremely sparse matrices. In particular, for an nxn matrix that has significantly fewer than n nonzeros these formats are rather wasteful because on top of the data and indices they carry a field indptr (delineating rows or columns) of size n+1.
Therefore assuming a properly utilized csc or csr matrix it is reasonable to expect row or column sums not to be sparse and the corresponding method should return a dense vector.
I'm aware that your question of "why" mostly targets the motivation behind the design decision, but anyway I tracked down how the result of csc_matrix.sum(axis=1) actually becomes a numpy matrix.
The csc_matrix class inherits from the _cs_matrix class which inherits from the _data_matrix class which inherits from the spmatrix base class. This last one implements .sum(ax) as
if axis == 0: # sum over columns ret = np.asmatrix(np.ones( (1, m), dtype=res_dtype)) * self else: # sum over rows ret = self * np.asmatrix(np.ones((n, 1), dtype=res_dtype)) In other words, as also noted in a comment, the column/row sums are computed by multiplying with a dense row or column matrix of ones, respectively. The result of this operation will be a dense matrix which you see on output.
While some of the subclasses override their .sum() method, as far as I could tell this only happens for the axis=None case, so the result which you see can be attributed to the above block of code.
The csr and csc formats were developed for linear algeba, especially the solution of large, but sparse, linear equations
A*x = b x = b/A A must be invertible, and can't have all 0's rows or columns.
A.sum(1) is done by matrix multiplication, with a (n,1) array of 1s.
With your mat:
In [203]: np.allclose(mat*np.mat(np.ones((120,1))), mat.sum(1)) Out[203]: True Doing that myself is actually a bit faster (overhead somewhere?)
In [204]: timeit mat.sum(1) 92.7 µs ± 111 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each) In [205]: timeit mat*np.mat(np.ones((120,1))) 59.2 µs ± 53.1 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each) I could also do this with a sparse matrix:
In [209]: mat*sparse.csc_matrix(np.ones((120,1))) Out[209]: <100x1 sparse matrix of type '<class 'numpy.float64'>' with 100 stored elements in Compressed Sparse Column format> In [211]: np.allclose(mat.sum(1),_.todense()) Out[211]: True But the time is slower, even if I move the sparse creation outside the loop:
In [213]: %%timeit I=sparse.csc_matrix(np.ones((120,1))) ...: mat*I ...: 215 µs ± 401 ns per loop (mean ± std. dev. of 7 runs, 1000 loops each) If mat was (115100,10) with lots of all 0 rows, this all sparse approach could give both time and space savings.
mat[:,:10] is also performed with matrix multiplication, with a sparse extractor matrix.
It is actually slower than the row sum:
In [247]: timeit mat[:,:10] 305 µs ± 10.4 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each) In [248]: timeit mat[:,:10].sum(1) 384 µs ± 9.05 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each) I can combine the column selection with sum using:
In [252]: I = sparse.lil_matrix((120,1),dtype=int); I[:10,:]=1; I=I.tocsc() In [253]: I Out[253]: <120x1 sparse matrix of type '<class 'numpy.int64'>' with 10 stored elements in Compressed Sparse Column format> In [254]: np.allclose((mat*I).todense(),mat[:,:10].sum(1)) Out[254]: True Timing on this mat*I is slower, though I could improve the I construction step.
I = sparse.csc_matrix((np.ones(10,int), np.arange(10), np.array([0,10])), shape=(120,1)) 
