How to invert vector of matrices in numpy?-CodePudding

Please assume a vector of invertible matrices:

import numpy as np


a = np.arange(120).reshape((2, 2, 5, 6))

I want to invert the matrices over their defined axes:

b = np.linalg.inv(a, axis1=0, axis2=1)

but this does not seems supported.

How to achieve this?

CodePudding user response：

If you know that the matrices are 2x2 you can do that easily using the standard formula for inverting such matrices; otherwise, I fear the only reasonable solution would be to do it with for loops? For example, the following works for any shape (modifying the sizes adequately):

b = np.stack([np.linalg.inv(a[:, :, i, j]) for i in range(a.shape[2]) for j in range(a.shape[3])], axis=2)
b = b.reshape(2, 2, 5, 6)

as checked by

for i in range(a.shape[2]):
    for j in range(a.shape[3]):
        assert np.allclose(np.dot(a[:,:,i,j], b[:,:,i,j]), np.eye(2))

In the specific 2x2 case you can do the following, which is fully vectorized hence probably faster:

determinants = a[0, 0] * a[1, 1] - a[0, 1] * a[1, 0]
b = 1 / determinants * np.stack([
    np.stack([a[1, 1], -a[0, 1]]),
    np.stack([-a[1, 0], a[0, 0]]),
])

On the specific (small) input size, the second solution is about 10 times faster in my tests (43us vs. 537us).

CodePudding user response：

inv docs specifies its array input as:

a : (..., M, M) array_like
    Matrix to be inverted.

You have a

a = np.arange(120).reshape((2, 2, 5, 6))
(M,M,...)

The dimensions are in the wrong order - change them!

In [44]: a = np.arange(120).reshape((2, 2, 5, 6))

Change the axes to the order that inv accepts:

In [45]: A = a.transpose(2,3,0,1)
In [46]: Ai = np.linalg.inv(A)
In [47]: Ai.shape
Out[47]: (5, 6, 2, 2)
In [48]: ai = Ai.transpose(2,3,0,1)    # and back
In [49]: ai.shape
Out[49]: (2, 2, 5, 6)

I was going to test the result, but got:

In [50]: x = a@ai
Traceback (most recent call last):
  File "<ipython-input-50-9dfe3616745d>", line 1, in <module>
    x = a@ai
ValueError: matmul: Input operand 1 has a mismatch in its core dimension 0, with gufunc signature (n?,k),(k,m?)->(n?,m?) (size 5 is different from 6)

Like inv, matmul treats the last 2 dimensions as the matrix, the first 2 as 'batch':

In [51]: x = A@Ai
In [52]: x[0,0]
Out[52]: 
array([[1., 0.],
       [0., 1.]])
In [53]: x[0,3]
Out[53]: 
array([[1.00000000e 00, 1.38777878e-17],
       [4.44089210e-16, 1.00000000e 00]])

We can do the equivalent with einsum:

In [55]: x = np.einsum('ijkl,jmkl->imkl',a,ai)
In [56]: x[:,:,0,0]
Out[56]: 
array([[1., 0.],
       [0., 1.]])

You might want to change the original specification to match the inv and matmul usage. It could make life easier for you. Also remember that in numpy the trailing dimensions are the inner most ones.