Fast inverse and transpose matrix in Python

Question

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?

Felixfeliza

Your A is not a matrix but a tensor. For a tensor it is not clear how to define an inverse or a transpose. If you want to inverse/transpose a 2-dim array of matrices you might want to look at numpy's tensorinv. – Motherwell

Answer 1

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)

Nomen

It works great! Do you have a similar solution in the case of 6x6 matrix? – Felixfeliza

Computing the inverse by means of an adjoint is only efficient for small matrices. Inverting 6x6 matrices (assuming no structure can be exploited) is much more expensive, and I can imagine the loop over n is no longer that much of a bottleneck in that scenario. – Nomen

Math question: is the cross product necessary for computing adjoints? – Whose

Answer 2

for the transpose:

testing a bit in ipython showed:

In [1]: import numpy
In [2]: x = numpy.ones((5,6,3,4))
In [3]: numpy.transpose(x,(0,1,3,2)).shape
Out[3]: (5, 6, 4, 3)

so you can just do

Atrans = numpy.transpose(A,(0,1,3,2))

to transpose the second and third dimensions (while leaving dimension 0 and 1 the same)

for the inversion:

the last example of http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.inv.html#numpy.linalg.inv

Inverses of several matrices can be computed at once:

from numpy.linalg import inv
a = np.array([[[1., 2.], [3., 4.]], [[1, 3], [3, 5]]])
>>> inv(a)
array([[[-2. ,  1. ],
        [ 1.5, -0.5]],
       [[-5. ,  2. ],
        [ 3. , -1. ]]])

So i guess in your case, the inversion can be done with just

Ainv = inv(A)

and it will know that the last two dimensions are the ones it is supposed to invert over, and that the first dimensions are just how you stacked your data. This should be much faster

speed difference

for the transpose: your method needs 3.77557015419 sec, and mine needs 2.86102294922e-06 sec (which is a speedup of over 1 million times)

for the inversion: i guess my numpy version is not high enough to try that numpy.linalg.inv trick with (n,n,3,3) shape, to see the speedup there (my version is 1.6.2, and the docs i based my solution on are for 1.8, but it should work on 1.8, if someone else can test that?)

Liquesce

Indeed this should work as of the latest numpy; but I don't know how fast it will be. This will simple call inv for each 3x3 matrix in the array; but the overhead of calling inv without exploiting its 3x3 nature of the matrices may be substantial; I suspect it would be slower than the solution I posted; but itd be nice to try. – Nomen

@EelcoHoogendoorn not certain, but could it be faster if you just hardcode a 3X3 inversion method like you do, but push it in C, and just scipy.weave.blitz on it? - Anyone with numpy 1.8 able to speedtest to compare my method with yours? – Liquesce

Answer 3

Numpy has the array.T properties which is a shortcut for transpose.

For inversions, you use np.linalg.inv(A).

Anhydrous

There is no array.I. That is a method on matrix object. For array you use np.linalg.inv(A). – Chigetai

Answer 4

As posted by wim A.I also works on matrix. e.g.

print (A.I)

Northing

Looks like this answer have similar point to your another answer here. You can edit your previous answer if you want to improve it. Or add new answer IF ONLY you have different point. Anyway, welcome to StackOverflow :) – Dight

Answer 5

for numpy-matrix object, use matrix.getI. e.g.

A=numpy.matrix('1 3;5 6')
print (A.getI())

Northing

Can you explain your answer? – Raffo