Complexity of matrix inversion in numpy

Question

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?

Answer 1

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

I'll assume you actually need to compute an inverse in your algorithm.¹ 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.

---

¹ 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.

Answer 2

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.