[SciPy-user] Invert a sparse matrix

[Dominique Orban]

On 10/4/07, Robin <robince at gmail.com> wrote:
> On 10/4/07, Alan G Isaac <aisaac at american.edu> wrote:
> >
> > Are you sure you need an inverse?
> > If you are just solving equation systems,
> > you usually do not.
>
> I'm pretty sure I need it... I am not solving equation systems. I am doing
> coordinate transformations in a large dimensional space (hence the large
> sparse matrix)... Some of the transformations required the inverse and the
> transpose of the core matrix.
>
> I'm actually trying to translate my code from MATLAB, where I used the
> sparse inverse, but I suppose if I want to do something like
> vec1 = A^(-T) ln (vec2)
> then I can use spsolve on A^T.vec1 = ln(vec2).
>
> The only thing is this is going to be inside an optimisation loop so I
> really need to get it as efficient as possible - I thought sparse matrix
> multiplication by the inverse is going to be a lot faster than solving each
> time...

I presume then that A does not change from one loop pass to the next.
It might be more efficient to factorize A once outside the loop and
keep its factors around (L and U if A is not symmetric, or just L if
it is symmetric and positive definite, etc). E.g., you can use UMFPACK
and call umfpack.solve, asking for a solve with A^T.

Cases where you *really* need the inverse are rare. Moreover, the
inverse of a sparse matrix is not necessarily sparse. Typically,
inverting is more expensive than an LU factorization and prone to
rounding errors.

"Almost anything you can do with A^{-1} can be done without it"
(George E. Forsythe and Cleve B. Moler, Computer Solution of Linear
Algebraic Systems, 1967).

See:
from scipy import linsolve
help(linsolve.umfpack)

Dominique