Inverting large sparse matrices with scipy

Question

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-diagonal elements (I'd rather not, but as a solution it'd be acceptable).

The matrices I need to invert are typically large(40000 *40000), and only have a handful of non-nonzero diagonals. My current approach is to build everything sparse, and then

posterior_covar = np.linalg.inv ( hessian.todense() )

this clearly takes a long time and plenty of memory.

Any hints, or it's just a matter of patience or making the problem smaller?

Answer 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]