Search code examples
pythonperformancenumpyscipysparse-matrix

Determining a sparse matrix quotient


I am looking to divide two sparce matricies in python 2.7, essentially k = numerator / denominator, with the result as a sparse matrix of type sp.csr_matrix. I am using scipy as sp and numpy as np.

To do this, I am following linear format of taking the dot product of numerator and inverse of denominator. Both items are of format sp.csr_matrix(([],([],[])),shape=[R,R]).

The calculation for k itself is

k = sp.csr_matrix(numerator.dot(sp.linalg.inv(denominator)))

Doing this returns the warning:

SparseEfficiencyWarning: splu requires CSC matrix format
warn('splu requires CSC matrix format', SparseEfficiencyWarning)

What does the above warning mean in relation to determining the identity of k as the quotient between two sparse matrices?

Is there a more efficient way to generate the dot product of a sparse matrix and a sparse matrix inverse in python (the quotient of two sparse matrices)?

I had previously found Inverting large sparse matrices with scipy, however I wonder if this may be outdated.


Solution

  • Borrowing from the 2013 answer:

    In [409]: a=np.random.rand(3,3)
    In [410]: A=sparse.csr_matrix(a)
    In [411]: np.linalg.inv(a)
    Out[411]: 
    array([[ 26.11275413,  -4.17749006,  -9.82626551],
           [-37.22611759,   9.38404027,  13.80073216],
           [  7.59314843,  -2.04314605,  -1.58410661]])
    

    The np inv is not sparse-aware:

    In [412]: np.linalg.inv(A) 
     ....
    LinAlgError: 0-dimensional array given. Array must be at least two-dimensional
    

    With from scipy.sparse import linalg:

    In [414]: linalg.inv(A).A
    /usr/lib/python3/dist-packages/scipy/sparse/linalg/dsolve/linsolve.py:243: SparseEfficiencyWarning: splu requires CSC matrix format
      warn('splu requires CSC matrix format', SparseEfficiencyWarning)
    /usr/lib/python3/dist-packages/scipy/sparse/linalg/dsolve/linsolve.py:161: SparseEfficiencyWarning: spsolve is more efficient when sparse b is in the CSC matrix format
      'is in the CSC matrix format', SparseEfficiencyWarning)
    Out[414]: 
    array([[ 26.11275413,  -4.17749006,  -9.82626551],
           [-37.22611759,   9.38404027,  13.80073216],
           [  7.59314843,  -2.04314605,  -1.58410661]])
    

    So use the csc format instead of csr:

    In [415]: A1=sparse.csc_matrix(a)
    In [416]: linalg.inv(A1).A
    Out[416]: 
    array([[ 26.11275413,  -4.17749006,  -9.82626551],
           [-37.22611759,   9.38404027,  13.80073216],
           [  7.59314843,  -2.04314605,  -1.58410661]])
    

    Same thing, but without the sparsity warning. Without getting into details, the inv must be using a method that iterates on columns rather than rows. It does spsolve(A, I) (I is a sparse eye matrix).