N*M*M张量的矢量化（部分）逆

Question

我和一年前问我的情况差不多： fast way to invert or dot kxnxn matrix

所以我有一个指数为a[n，I，j]的张量，它的维数为（n，M，M），我想把n中每个n的M*M平方矩阵部分求逆

例如，假设我有

In [1]:    a = np.arange(12)
           a.shape = (3,2,2)
           a

Out[1]: array([[[ 0,  1],
                  [ 2,  3]],

                  [[ 4,  5],
                  [ 6,  7]],

                  [[ 8,  9],
                  [10, 11]]])

那么for循环反转将如下所示：

^{pr2}$

根据github上的this issue，这显然将在numpy2.0中实现。。。在

我想我需要像seberg在github问题线程中提到的那样安装dev版本，但是现在有没有其他方法可以以矢量化的方式安装dev版本？在

Answer 1

更新： 在nump1.8和更高版本中，numpy.linalg中的函数是广义通用函数。 也就是说你现在可以这样做了：

import numpy as np
a = np.random.rand(12, 3, 3)
np.linalg.inv(a)

这将反转每个3x3数组，并以12x3x3数组的形式返回结果。 参见numpy 1.8 release notes。在

---

原始答案：

由于N相对较小，我们一次手动计算所有矩阵的LU分解如何。 这样可以确保所涉及的for循环相对较短。在

下面是如何使用普通NumPy语法完成此操作：

^{pr2}$

正如我所写的，pylu3d就地修改A来计算LU分解。 将每个NxN矩阵替换为其LU分解后，pylusolve可用于求解代表矩阵系统右侧的MxN数组b。 它对b进行适当的修改，并进行适当的反向替换以解决系统问题。 正如它所写的，这个实现不包括旋转，所以它在数值上不稳定，但在大多数情况下它应该能很好地工作。在

根据数组在内存中的排列方式，使用Cython可能还是快一点。 这里有两个Cython函数执行相同的操作，但是它们首先沿着M进行迭代。 它没有矢量化，但速度相对较快。在

from numpy cimport ndarray as ar
cimport cython

@cython.boundscheck(False)
@cython.wraparound(False)
def lu3d(ar[double,ndim=3] A):
    cdef int n, i, j, k, N=A.shape[0], h=A.shape[1], w=A.shape[2]
    for n in xrange(N):
        for j in xrange(h-1):
            for i in xrange(j+1,h):
                #change to L
                A[n,i,j] /= A[n,j,j]
                #change to U
                for k in xrange(j+1,w):
                    A[n,i,k] -= A[n,i,j] * A[n,j,k]

@cython.boundscheck(False)
@cython.wraparound(False)
def lusolve(ar[double,ndim=3] A, ar[double,ndim=2] b):
    cdef int n, i, j, N=A.shape[0], h=A.shape[1]
    for n in xrange(N):
        for j in xrange(h-1):
            for i in xrange(j+1,h):
                b[n,i] -= A[n,i,j] * b[n,j]
        for j in xrange(h-1,-1,-1):
            b[n,j] /= A[n,j,j]
            for i in xrange(j):
                b[n,i] -= A[n,i,j] * b[n,j]

你也可以尝试使用Numba，虽然在这个例子中我不能让它像Cython一样快。在