ale*_*lex 5 numpy matrix scipy
有没有一种使用numpy计算kxnxn矩阵的逆的快速方法(在每个k切片计算的倒数)？换句话说,有没有办法矢量化以下代码:
>>>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,:,:])
首先是关于求逆。我已经研究过np.linalg.tensorinv和np.linalg.tensorsolve。
我认为不幸的tensorinv是不会给你你想要的。它需要数组是“方形”的。这排除了你想要做的事情，因为他们对正方形的定义是，np.prod(a[:i]) == np.prod(a[i:])其中i是 0、1 或 2（一般是数组的轴之一）；这可以作为ind的第三个参数给出tensorinv。这意味着，如果您有一个长度为 M 的 NxN 矩阵的通用数组，则需要有例如（对于 i = 1）NxN == NxM，这通常是不正确的（在您的示例中实际上是正确的，但它确实无论如何都没有给出正确答案）。
现在，也许有些事情是可能的tensorsolve。a然而，在将矩阵数组作为第一个参数传递给 之前，这将涉及对矩阵数组进行一些繁重的构造工作tensorsolve。因为我们希望b成为“矩阵数组方程” a*b = 1（其中1是单位矩阵数组）的解，并且1具有与a和相同的形状b，所以我们不能简单地将a上面定义的 you 提供为 的第一个参数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.einsumbetter。例如：
a=np.random.random((7,2,2))
b=np.random.random((7,2,2))
np.einsum('ijk,ikl->ijl', a,b)
a这将对数组和中的所有 7 个“矩阵”进行矩阵乘法b。它似乎比上面博客文章中的数组方法慢大约 2 倍，但它仍然比您的示例中使用 for 循环快大约 70 倍。事实上，对于更大的数组（例如 10000 个 5x5 矩阵），该einsum方法似乎稍微快一些（不知道为什么）。
希望有帮助。