Python中的快速逆和转置矩阵

Question

我有一个大的矩阵A状的(n, n, 3, 3)与n约5000.现在我想找到矩阵的逆和转置A:

import numpy as np
A = np.random.rand(1000, 1000, 3, 3)
identity = np.identity(3, dtype=A.dtype)
Ainv = np.zeros_like(A)
Atrans = np.zeros_like(A)
for i in range(1000):
    for j in range(1000):
        Ainv[i, j] = np.linalg.solve(A[i, j], identity)
        Atrans[i, j] = np.transpose(A[i, j])

有没有更快,更有效的方法来做到这一点？

Answer 1

这是从我的一个项目中获得的,我也在许多3x3矩阵上进行矢量化线性代数.

请注意,只有一个超过3的循环; 不是n上的循环,因此代码在重要维度中被矢量化.我不想保证这与C/numba扩展相比如何做同样的事情,性能明智.这可能会大大加快,但至少这会使循环超过水中的n.

def adjoint(A):
    """compute inverse without division by det; ...xv3xc3 input, or array of matrices assumed"""
    AI = np.empty_like(A)
    for i in xrange(3):
        AI[...,i,:] = np.cross(A[...,i-2,:], A[...,i-1,:])
    return AI

def inverse_transpose(A):
    """
    efficiently compute the inverse-transpose for stack of 3x3 matrices
    """
    I = adjoint(A)
    det = dot(I, A).mean(axis=-1)
    return I / det[...,None,None]

def inverse(A):
    """inverse of a stack of 3x3 matrices"""
    return np.swapaxes( inverse_transpose(A), -1,-2)
def dot(A, B):
    """dot arrays of vecs; contract over last indices"""
    return np.einsum('...i,...i->...', A, B)


A = np.random.rand(2,2,3,3)
I = inverse(A)
print np.einsum('...ij,...jk',A,I)

Answer 2

对于转置：

在 ipython 中进行一些测试表明：

In [1]: import numpy
In [2]: x = numpy.ones((5,6,3,4))
In [3]: numpy.transpose(x,(0,1,3,2)).shape
Out[3]: (5, 6, 4, 3)

所以你可以这样做

Atrans = numpy.transpose(A,(0,1,3,2))

转置第二个和第三个维度（同时保持维度 0 和 1 相同）

对于反转：

http://docs.scipy.org/doc/numpy/reference/ generated/numpy.linalg.inv.html#numpy.linalg.inv的最后一个例子

可以同时计算多个矩阵的逆：

from numpy.linalg import inv
a = np.array([[[1., 2.], [3., 4.]], [[1, 3], [3, 5]]])
>>> inv(a)
array([[[-2. ,  1. ],
        [ 1.5, -0.5]],
       [[-5. ,  2. ],
        [ 3. , -1. ]]])

所以我想在你的情况下，反转可以通过

Ainv = inv(A)

它会知道最后两个维度是它应该反转的维度，而第一个维度正是您堆叠数据的方式。这应该会快得多

速度差

对于转置：你的方法需要 3.77557015419 秒，而我的方法需要 2.86102294922e-06 秒（加速超过 100 万倍）

对于反转：我想我的 numpy 版本不够高，无法尝试使用 (n,n,3,3) 形状的 numpy.linalg.inv 技巧，以查看那里的加速（我的版本是 1.6.2，文档我的解决方案基于 1.8，但它应该适用于 1.8，如果其他人可以测试的话？）

Answer 3

Numpy 具有转置array.T快捷方式的属性。

对于反转，您可以使用np.linalg.inv(A).