Search code examples
pythonnumpymatrixsparse-matrixdiagonal

numpy functions that can/cannot be used on a compressed sparse row (CSR) matrix


I am a newbie in Python and I have a (probably very naive) question. I have a CSR (compressed sparse row) matrix to work on (let's name it M), and looks like some functions that are designed for a 2d numpy array manipulation work for my matrix while some others do not.

For example, numpy.sum(M, axis=0) works fine while numpy.diagonal(M) gives an error saying {ValueError}diag requires an array of at least two dimensions.

So is there a rationale behind why one matrix function works on M while the other does not?

And a bonus question is, how to get the diagonal elements from a CSR matrix given the above numpy.diagonal does not work for it?


Solution

  • The code for np.diagonal is:

    return asanyarray(a).diagonal(offset=offset, axis1=axis1, axis2=axis2)
    

    That is, it first tries to turn the argument into an array, for example if it is a list of lists. But that isn't the right way to turn a sparse matrix into a ndarray.

    In [33]: from scipy import sparse                                               
    In [34]: M = sparse.csr_matrix(np.eye(3))                                       
    In [35]: M                                                                      
    Out[35]: 
    <3x3 sparse matrix of type '<class 'numpy.float64'>'
        with 3 stored elements in Compressed Sparse Row format>
    In [36]: M.A                                  # right                                  
    Out[36]: 
    array([[1., 0., 0.],
           [0., 1., 0.],
           [0., 0., 1.]])
    In [37]: np.asanyarray(M)                    # wrong                           
    Out[37]: 
    array(<3x3 sparse matrix of type '<class 'numpy.float64'>'
        with 3 stored elements in Compressed Sparse Row format>, dtype=object)
    

    The correct way to use np.diagonal is:

    In [38]: np.diagonal(M.A)                                                       
    Out[38]: array([1., 1., 1.])
    

    But no need for that. M already has a diagonal method:

    In [39]: M.diagonal()                                                           
    Out[39]: array([1., 1., 1.])
    

    np.sum does work, because it delegates the action to a method (look at its code):

    In [40]: M.sum(axis=0)                                                          
    Out[40]: matrix([[1., 1., 1.]])
    In [41]: np.sum(M, axis=0)                                                      
    Out[41]: matrix([[1., 1., 1.]])
    

    As a general rule, try to use sparse functions and methods on sparse matrices. Don't count on numpy functions to work right. sparse is built on numpy, but numpy does not 'know' about sparse.