How to calculate the generalized inverse of a Sparse Matrix in scipy

Question

I have a sparse matrix W, when I use linalg.pinv(W), it throws some errors:

Traceback (most recent call last):
  File "/Users/ad9075/PycharmProjects/bednmf/test.py", line 14, in testNmfRun
    self.factor = factorization(self.V)
  File "/Users/ad9075/PycharmProjects/bednmf/nmf.py", line 18, in factorization
    W_trans = linalg.pinv(W)
  File "/Library/Python/2.7/site-packages/scipy/linalg/basic.py", line 540, in pinv
    b = np.identity(a.shape[0], dtype=a.dtype)
IndexError: tuple index out of range`

But when I modify it to linalg.pinv(W.todense()), it works well. However, do I really need to convert the sparse matrix if I want to calculate the generaized inverse? Does anyone have ideas about this?

Thanks!

Answer 1

The inverse (and generalized inverse) of a sparse matrix is usually dense, unless you can permute the rows and columns of the matrix so that it becomes block diagonal.

So your problem splits into two parts: (i) find a permutation that makes it block-diagonal, and (ii) compute the generalized inverse using linalg.pinv separately for each block. If your matrix is small enough, just converting it to a dense matrix first and then computing the pseudoinverse is also efficient.

If you on the other hand want to compute something like "A^{-1} x", using gmres or some other iterative routine may be a more efficient solution.

Answer 2

Your answer lies here : Good matrix inversion routines in C ( don't let the title disturb you, the discussion is about scipy )

But you have to know the inverse matrix is not necessarily a sparse matrix and, as said by George E. Forsythe and Cleve B. Moler :

"Almost anything you can do with A^{-1} can be done without it"