CODEWITHSUNDEEP
pythonnumpylinear-algebra
Freddie Witherden
Freddie Witherden

Reputation: 2426

Efficient 3x3 and 2x2 Determinants in NumPy

I have a large numpy array, arr, of shape (N, D, M, D) where D is either two or three. The array can be thought of as a block of (D,D) matrices blocked together in the N and M dimensions. I wish to take the determinant of each of these (D,D) matrices:

out = np.empty((N,M))
for i in xrange(N):
    for j in xrange(M):
        out[i,j] = np.linalg.det(arr[i,:,j,:])

The only problem is that np.linalg.det does not special-case small matrices and so does a complete BLAS call and LU factorization for each block. Indeed, for a 2 by 2 matrix np.linalg.det takes ~40us on my recent Core i7 system. What are my options for improving the performance of this fragment?

Upvotes: 1

Views: 1280

Answers (2)

senderle
senderle

Reputation: 151017

To clarify duffymo's point, you can easily exploit the block structure with your own function, using vectorized numpy operations. Whether this gives you a speedup may depend on how large M and N are:

>>> a = numpy.arange(2 * 2 * 2 * 2).reshape(2, 2, 2, 2)
>>> a[:,0,:,0] * a[:,1,:,1] - a[:,1,:,0] * a[:,0,:,1]
array([[-4, -4],
       [-4, -4]])
>>> a = numpy.arange(2 * 3 * 2 * 3).reshape(2, 3, 2, 3)
>>> a[:,0,:,0] * a[:,1,:,1] * a[:,2,:,2] + \
    a[:,0,:,1] * a[:,1,:,2] * a[:,2,:,0] + \
    a[:,0,:,2] * a[:,1,:,0] * a[:,2,:,1] - \
    a[:,0,:,0] * a[:,1,:,2] * a[:,2,:,1] - \
    a[:,0,:,1] * a[:,1,:,0] * a[:,2,:,2] - \
    a[:,0,:,2] * a[:,1,:,1] * a[:,2,:,0]
array([[0, 0],
       [0, 0]])

Upvotes: 3

duffymo
duffymo

Reputation: 308823

You don't need LU decomposition to calculate the determinant of 2x2 or 3x3 matricies. The formulas are well known. Why not write your own functions and call them?

Upvotes: 1

Related Questions