2

As part of a python code for a numerical calculation, I must invert somewhat large (sparse) matrices (~100x100) many times. I would really like to speed up the program, and one of the ways suggested to me is to call to a subroutine in C for the matrix inversion step.

Are there any recommended efficient and well-tested C routines for this task?

Thank you.

Shai
  • 111,146
  • 38
  • 238
  • 371
AKC
  • 21
  • 2
  • 2
    Have you looked at the NumPy + SciPy modules? – adelbertc Jun 25 '12 at 05:51
  • Yes, NumPy is actually what I am using right now, and each data point takes around 5s to do, which is too slow considering that I will need lots of points. I've read that NumPy is considered fairly efficient already, but would calling something from a compiled language like C save me a lot of time? Thank you for your comment! – AKC Jun 25 '12 at 05:55
  • 1
    NumPy and SciPy both have much of their time-critical loops written in C and/or FORTRAN - not having much experience with making Python interface with C myself, I'm not sure how much improvement you can get from writing one yourself - matrix inversion tends to be an expensive computation. That being said, a popular matrix inversion algorithm seems to be Gaussian Elimination. – adelbertc Jun 25 '12 at 06:13
  • Oh I see now. I guess I will have to settle with this speed for now and be more clever in sampling only the more important domains. Still, it is good to learn that NumPy routines are already one of the more optimized choices - Python is just so easy to learn and use! Thanks again! – AKC Jun 25 '12 at 06:25
  • You can check out [CSPARSE](http://people.sc.fsu.edu/~jburkardt/c_src/csparse/csparse.html)... – Eitan T Jun 25 '12 at 07:51
  • Gauss-Jordan method is a pretty popular algorithm that inverses a matrix. it's really easy to implement. I don't know if you'll save any time though, but when I get some free time, I'll implement it. (cause it's also fun!) – moldovean Mar 23 '15 at 14:38

1 Answers1

2
>>> from numpy import *
>>> from numpy.linalg import inv
>>> from scipy.sparse import csr_matrix
>>> m = matrix([[3,1,5],[1,0,8],[2,1,4]])
>>> s = csr_matrix(m)
>>> invs = inv(a) # Inverse sparse matrix
>>> dot(a,inva) # Check the result, should be eye(3) within machine precision
csr_matrix([[ 1.00000000e-00, 2.77555756e-17, 3.60822483e-16],
           [ 0.00000000e+00, 1.00000000e+00, 0.00000000e+00],
           [ -1.11022302e-16, 0.00000000e+00, 1.00000000e+00]])

Is it really the inverse you need? You may be able to achieve your goal without inversion:

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.

-- http://mail.scipy.org/pipermail/scipy-user/2007-October/013936.html

--> http://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.factorized.html#scipy.sparse.linalg.factorized

Thomas
  • 11,757
  • 4
  • 41
  • 57