Is there a way to effectively invert an array of matrices with numpy?

Question

I usually inverted an array of 3x3 matrices in a for loop, as in the example below. Unfortunately, for loops are slow. Is there a faster and more efficient way to do this?

 import numpy as np A = np.random.rand(3,3,100) Ainv = np.zeros_like(A) for i in range(100): Ainv[:,:,i] = np.linalg.inv(A[:,:,i]) 

Answer 1

Turns out you burn two levels in numpy.linalg code. If you look at numpy.linalg.inv, you can see it simply as a call to numpy.linalg.solve (A, inv (A.shape [0]). This leads to the reconstruction of the identity matrix at each iteration of your loop. Since all your arrays are the same size, it’s a waste of time. Skipping this step by pre-allocating the identification matrix will hit ~ 20% of the time (fast_inverse). My testing assumes that pre-allocating the array or selecting it from the list of results does not really matter.

Look one level deeper and you will find the challenge of a routine program, but it wraps up a few sanity checks. If you cross out all of this and just call lapack in your for loop (since you already know the size of your matrix and maybe you know that it is real and not complicated), everything works much faster (note that I made my array more) :

 import numpy as np A = np.random.rand(1000,3,3) def slow_inverse(A): Ainv = np.zeros_like(A) for i in range(A.shape[0]): Ainv[i] = np.linalg.inv(A[i]) return Ainv def fast_inverse(A): identity = np.identity(A.shape[2], dtype=A.dtype) Ainv = np.zeros_like(A) for i in range(A.shape[0]): Ainv[i] = np.linalg.solve(A[i], identity) return Ainv def fast_inverse2(A): identity = np.identity(A.shape[2], dtype=A.dtype) return array([np.linalg.solve(x, identity) for x in A]) from numpy.linalg import lapack_lite lapack_routine = lapack_lite.dgesv # Looking one step deeper, we see that solve performs many sanity checks. # Stripping these, we have: def faster_inverse(A): b = np.identity(A.shape[2], dtype=A.dtype) n_eq = A.shape[1] n_rhs = A.shape[2] pivots = zeros(n_eq, np.intc) identity = np.eye(n_eq) def lapack_inverse(a): b = np.copy(identity) pivots = zeros(n_eq, np.intc) results = lapack_lite.dgesv(n_eq, n_rhs, a, n_eq, pivots, b, n_eq, 0) if results['info'] > 0: raise LinAlgError('Singular matrix') return b return array([lapack_inverse(a) for a in A]) %timeit -n 20 aI11 = slow_inverse(A) %timeit -n 20 aI12 = fast_inverse(A) %timeit -n 20 aI13 = fast_inverse2(A) %timeit -n 20 aI14 = faster_inverse(A) 

The results are impressive:

 20 loops, best of 3: 45.1 ms per loop 20 loops, best of 3: 38.1 ms per loop 20 loops, best of 3: 38.9 ms per loop 20 loops, best of 3: 13.8 ms per loop 

EDIT: I did not look closely enough at what was returning in the decision process. It turns out that the matrix 'b' is overwritten and contains the result at the end. This code now gives consistent results.

Answer 2

For loops, it really is not necessarily much slower than the alternatives, in which case this will not help you. But here is the suggestion:

 import numpy as np A = np.random.rand(100,3,3) #this is to makes it #possible to index #the matrices as A[i] Ainv = np.array(map(np.linalg.inv, A)) 

The timing of this solution compared to your solution gives a small but noticeable difference:

 # The for loop: 100 loops, best of 3: 6.38 ms per loop # The map: 100 loops, best of 3: 5.81 ms per loop 

I tried to use the vectorize numpy procedure with the hope of creating an even cleaner solution, but I have to figure it out. The reordering in array A is probably the most significant change since it exploits the fact that numpy arrays are sorted by columns, and therefore linear data reading is so fast.