python - Python 中的快速逆矩阵和转置矩阵

Question

我有一个形状为 (n, n, 3, 3) 的大型矩阵 A，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

这取 self 的一个项目，我也在该项目中对许多 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)