快速反转或点积kxnxn矩阵的方法

Question

是否有一种快速的方法，使用numpy计算kxnxn矩阵的逆（在每个k-slice上计算逆）？换句话说，是否有一种方法可以将以下代码进行向量化：

>>>from numpy.linalg import inv
>>>a-random(4*2*2).reshape(4,2,2)
>>>b=a.copy()
>>>for k in range(len(a)):
>>>    b[k,:,:] = inv(a[k,:,:])

- alex

---

我觉得 numpy.linalg.tensorinv 应该可以做到这一点，但我还没完全理解它…… - Joe Kington

Answer 1

首先是求逆的问题。我研究了 np.linalg.tensorinv 和 np.linalg.tensorsolve 两种方法。

不幸的是，tensorinv 方法不能给你想要的结果。它需要数组是“方阵”的。这将排除你想做的事情，因为他们定义的方阵是 i 是 0、1 或 2（一般情况下数组的一个轴）时，np.prod(a[:i]) == np.prod(a[i:])，可以作为 tensorinv 的第三个参数 ind 给出。这意味着如果你有一个长度为 M 的 NxN 矩阵的通用数组，则需要例如（对于 i = 1）NxN == NxM，这在一般情况下不成立（在你的示例中确实成立，但无论如何都不能得到正确答案）。

现在，也许使用 tensorsolve 可能会有所帮助。然而，这将涉及到在将 a 矩阵数组作为 tensorsolve 的第一个参数传递之前对其进行一些重建工作。因为我们希望 b 是“矩阵数组方程”的解 a*b = 1 （其中 1 是恒等矩阵数组）且 1 拥有与 a 和 b 相同的形状，我们不能简单地将你上面定义的 a 作为 tensorsolve 的第一个参数提供。相反，它需要是形状为（M，N，N，M，N，N）、（M，N，N，N，M，N）或（M，N，N，N，N，M）的数组。这是必要的，因为 tensorsolve 将在最后三个轴上乘以 b 并对它们求和，以使结果（该函数的第二个参数）再次具有形状（M，N，N）。

其次是点积问题（标题似乎也是你问题的一部分）。这很可行。有两种方法。

首先，James Hensman 的这篇博客文章 给出了一些不错的建议。

其次，我个人更喜欢使用 np.einsum，因为更清晰易懂。例如：

a=np.random.random((7,2,2))
b=np.random.random((7,2,2))
np.einsum('ijk,ikl->ijl', a,b)

这将会对数组a和b中的7个“矩阵”进行矩阵乘法运算。它似乎比上面博客文章中介绍的数组方法慢2倍左右，但仍然比使用for循环（就像你的示例）快70倍左右。实际上，对于更大的数组（例如10000个5x5矩阵），einsum方法似乎略快（不确定为什么）。希望这能帮到你。

- egpbos