CODEWITHSUNDEEP
pythonnumpytransposematrix-inverse
user2863620
user2863620

Reputation: 643

Fast inverse and transpose matrix in Python

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?

Upvotes: 6

Views: 27035

Answers (5)

Sunil Patil
Sunil Patil

Reputation: 21

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

print (A.I)

Upvotes: 1

Sunil Patil
Sunil Patil

Reputation: 21

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

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

Upvotes: 0

usethedeathstar
usethedeathstar

Reputation: 2309

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?)

Upvotes: 5

Eelco Hoogendoorn
Eelco Hoogendoorn

Reputation: 10759

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)

Upvotes: 7

Louis Thibault
Louis Thibault

Reputation: 21410

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

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

Upvotes: 3

Related Questions