如何在numpy中反转矩阵的体积？

Question

请假设一个可逆矩阵向量：

import numpy as np


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

我想在它们定义的轴上反转矩阵：

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

但这似乎不受支持。

如何做到这一点？

Answer 1

inv docs 将其数组输入指定为：

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

你有一个

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

尺寸顺序错误 - 更改它们！

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

将坐标轴更改为inv 接受的顺序：

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)

我正要测试结果，但得到了：

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)

与inv 一样，matmul 将最后 2 个维度视为 matrix，将前 2 个维度视为“批次”：

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

我们可以用einsum做同样的事情：

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

您可能希望更改原始规范以匹配 inv 和 matmul 用法。它可以让你的生活更轻松。还要记住，在numpy 中，尾随维度是最里面的维度。

Answer 2

如果您知道矩阵是 2x2，那么您可以使用标准公式轻松地对此类矩阵求逆；否则，我担心唯一合理的解决方案是使用 for 循环？例如，以下适用于任何形状（充分修改大小）：

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)

经检查

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

在特定的 2x2 情况下，您可以执行以下操作，这是完全矢量化的，因此可能更快：

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]]),
])

在特定（小）输入大小上，第二种解决方案在我的测试中大约快 10 倍（43us 与 537us）。