Inverting Large Sparse Matrices With Scipy

I have to invert a large sparse matrix. I cannot escape from the matrix inversion, the only shortcut would be to just get an idea of the main diagonal elements, and ignore the off-

Solution 1:

I don't think that the sparse module has an explicit inverse method, but it does have sparse solvers. Something like this toy example works:

>>> a = np.random.rand(3, 3)
>>> a
array([[ 0.31837307,  0.11282832,  0.70878689],
       [ 0.32481098,  0.94713997,  0.5034967 ],
       [ 0.391264  ,  0.58149983,  0.34353628]])
>>> np.linalg.inv(a)
array([[-0.29964242, -3.43275347,  5.64936743],
       [-0.78524966,  1.54400931, -0.64281108],
       [ 1.67045482,  1.29614174, -2.43525829]])

>>> a_sps = scipy.sparse.csc_matrix(a)
>>> lu_obj = scipy.sparse.linalg.splu(a_sps)
>>> lu_obj.solve(np.eye(3))
array([[-0.29964242, -0.78524966,  1.67045482],
       [-3.43275347,  1.54400931,  1.29614174],
       [ 5.64936743, -0.64281108, -2.43525829]])

Note that the result is transposed!

If you expect your inverse to also be sparse, and the dense return from the last solve won't fit in memory, you can also generate it one row (column) at a time, extract the non-zero values, and build the sparse inverse matrix from those:

>>>for k in xrange(3) :...    b = np.zeros((3,))...    b[k] = 1...print lu_obj.solve(b)... 
[-0.29964242 -0.78524966  1.67045482]
[-3.43275347  1.54400931  1.29614174]
[ 5.64936743 -0.64281108 -2.43525829]