I have a block-diagonal sparse matrix in python with hundreds of thousands to millions of rows, but with many different sized blocks between 1x1 and 6x6; for example plt.spy() on a submatrix:

I need the inverse of this matrix without looping in python (too slow). I should be able to extract the block diagonal to get a list of 2d arrays (have not actually implemented this yet), but even then I cannot figure out how to apply, e.g., np.linalg.inv() to a bunch of different size matrices. Any ideas on inverting in sparse matrix form or via a list of diagonal blocks??

CodePudding user response：

AFAIK, there is no way to do that efficiently using Numpy. The best solution would be to group the blocks by size to compute them in a vectorized way (most Numpy function cannot work on array of different sizes). But this solution is not great because Numpy is not designed to compute small arrays and the overhead will be pretty big for such a small blocks.

A solution is to use Numba so to generate a fast native code to compute your blocks efficiently. Numba supports typed lists and np.linalg.inv. If your input is a sparse matrix the best solution would be to extract the sub-blocks directly from it. However, Numba does not supports the Scipy's sparse matrices yet. That being said, you can for example extract the data and indices fields (which are Numpy array) of a CSR matrix so to compute them using Numba although this is a bit more complex to do.