Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I am solving differential equations that require to invert dense square matrices. This matrix inversion consumes the most of my computation time, so I was wondering if I am using the fastest algorithm available.

My current choice is numpy.linalg.inv. From my numerics I see that it scales as $O(n^3)$ where n is the number of rows, so the method seems to be Gaussian elimination.

According to Wikipedia, there are faster algorithms avaliable. Does anyone know if there is a library that implements these?

I wonder, why isn't numpy using these faster algorithms?

share|improve this question
    
You need to perform your matrices before. Look at Scipy. Sparse for your help. It contains many tools you need. – Tobal Feb 11 at 18:34
    
@Tobal not sure I follow... how would you "perform" a matrix? and exactly how would scipy.sparse help? – GoHokies Feb 11 at 19:24
    
@GoHokies scipy is a complement to numpy. Dense/sparse matrices must be implemented well before you do some calculations, it improves your calculations. Please read this docs.scipy.org/doc/scipy/reference/sparse.html it explains best than me with my bad english. – Tobal Feb 11 at 20:00
    
@Tobal The question specifically refers to dense matrices, so I don't see how scipy.sparse is relevant here? – Christian Clason Feb 11 at 20:21
2  
@Tobal -- I think I still don't understand. What exactly do you mean with "preform your matrices", and "matrices must be implemented well before you do some calculations"? Regarding your last comment, surely you will agree that the techniques that can be used for sparse and dense matrices are very different. – Wolfgang Bangerth Feb 11 at 23:23
up vote 14 down vote accepted

(This is getting too long for comments...)

I'll assume you actually need to compute an inverse in your algorithm.1 First, it is important to note that these alternative algorithms are not actually claimed to be faster, just that they have better asymptotic complexity (meaning the required number of elementary operations grows more slowly). In fact, in practice these are actually (much) slower than the standard approach (for given $n$), for the following reasons:

  1. The $\mathcal{O}$-notation hides a constant in front of the power of $n$, which can be astronomically large -- so large that $C_1 n^3$ can be much smaller than $C_2 n^{2.x}$ for any $n$ that can be handled by any computer in the foreseeable future. (This is the case for the Coppersmith–Winograd algorithm, for example.)

  2. The complexity assumes that every (arithmetical) operation takes the same time -- but this is far from true in actual practice: Multiplying a bunch of numbers with the same number is much faster than multiplying the same amount of different numbers. This is due to the fact that the major bottle-neck in current computing is getting the data into cache, not the actual arithmetical operations on that data. So an algorithm which can be rearranged to have the first situation (called cache-aware) will be much faster than one where this is not possible. (This is the case for the Strassen algorithm, for example.)

Also, numerical stability is at least as important as performance; and here, again, the standard approach usually wins.

For this reason, the standard high-performance libraries (BLAS/LAPACK, which Numpy calls when you ask it to compute an inverse) usually only implement this approach. Of course, there are Numpy implementations of, e.g., Strassen's algorithm out there, but an $\mathcal{O}(n^3)$ algorithm hand-tuned at assembly level will soundly beat an $\mathcal{O}(n^{2.x})$ algorithm written in a high-level language for any reasonable matrix size.


1 But I'd be amiss if I didn't point out that this is very rarely really necessary: anytime you need to compute a product $A^{-1}b$, you should instead solve the linear system $Ax=b$ (e.g., using numpy.linalg.solve) and use $x$ instead -- this is much more stable, and can be done (depending on the structure of the matrix $A$) much faster. If you need to use $A^{-1}$ multiple times, you can precompute a factorization of $A$ (which is usually the most expensive part of the solve) and reuse that later.

share|improve this answer
    
Great answer, thank you sir, in particular for pointing out the devil in the details (constants in big O notation) that makes a big difference between theoretical speed and practical speed. – gaborous Feb 11 at 20:23
    
I think that the "inverse is rarely necessary" part should be emphasized more. If the purpose is to solve a system of differential equations, it does not seem likely that a full inverse is needed. – Jared Goguen Feb 12 at 2:15
    
@o_o Well, that was my first original comment (which I deleted after consolidating them all into one answer). But I thought, for the benefit of the site (and later readers), an answer should answer the actual question in the question (which is both reasonable and on-topic), even if there's an XY problem behind it. Also, I didn't want to sound too admonishing... – Christian Clason Feb 12 at 8:24
    
Thank you, I was not aware that the constats were that high. That explains it well and I can stick to numpy without thinking it could be done much faster. – physicsGuy Feb 12 at 9:30
1  
As I wrote, you can in almost all cases rewrite your algorithm to replace operations involving the inverse with solving the corresponding linear system (or in this case, sequence of linear systems) -- if you are interested, you could ask a separate question about that ("can I avoid inverting matrices in this algorithm?"). And yes, since the number of matrices does not depend on $n$, the complexity is still the same (you just get a bigger constant -- by a factor of four in your case). – Christian Clason Feb 12 at 10:08

You should probably note that, buried deep inside the numpy source code (see https://github.com/numpy/numpy/blob/master/numpy/linalg/umath_linalg.c.src) the inv routine attempts to call the dgetrf function from your system LAPACK package, which then performs an LU decomposition of your original matrix. This is morally equivalent to Gaussian elimination, but can be tuned to a slightly lower complexity by using faster matrix multiplication algorithms in a high-performance BLAS.

If you follow this route, you should be warned that forcing the entire library chain to use the new library rather than the system one which came with your distribution is fairly complex. One alternative on modern computer systems is to look at parallelized methods using packages like scaLAPACK or (in the python world) petsc4py. However these are typically happier being used as iterative solvers for linear algebra systems than applied to direct methods and PETSc in particular targets sparse systems more than dense ones.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.