Fast inverse and transpose matrix in Python

matrix-inverse, numpy, python, transpose

Solution

This is taken from a project of mine, where I also do vectorized linear algebra on many 3x3 matrices.

Note that there is only a loop over 3; not a loop over n, so the code is vectorized in the important dimensions. I don't want to vouch for how this compares to a C/numba extension to do the same thing though, performance wise. This is likely to be substantially faster still, but at least this blows the loops over n out of the water.

def adjoint(A):
    """compute inverse without division by det; ...xv3xc3 input, or array of matrices assumed"""
    AI = np.empty_like(A)
    for i in xrange(3):
        AI[...,i,:] = np.cross(A[...,i-2,:], A[...,i-1,:])
    return AI

def inverse_transpose(A):
    """
    efficiently compute the inverse-transpose for stack of 3x3 matrices
    """
    I = adjoint(A)
    det = dot(I, A).mean(axis=-1)
    return I / det[...,None,None]

def inverse(A):
    """inverse of a stack of 3x3 matrices"""
    return np.swapaxes( inverse_transpose(A), -1,-2)
def dot(A, B):
    """dot arrays of vecs; contract over last indices"""
    return np.einsum('...i,...i->...', A, B)


A = np.random.rand(2,2,3,3)
I = inverse(A)
print np.einsum('...ij,...jk',A,I)

Problem

I have a large matrix `A` of shape `(n, n, 3, 3)` with `n` is about `5000`. Now I want find the inverse and transpose of matrix `A`: ``` import numpy as np A = np.random.rand(1000, 1000, 3, 3) identity = np.identity(3, dtype=A.dtype) Ainv = np.zeros_like(A) Atrans = np.zeros_like(A) for i in range(1000): for j in range(1000): Ainv[i, j] = np.linalg.solve(A[i, j], identity) Atrans[i, j] = np.transpose(A[i, j]) ``` Is there a faster, more efficient way to do this?

Original source