NumPy 中的高效 3x3 和 2x2 行列式答案

【问题标题】：Efficient 3x3 and 2x2 Determinants in NumPyNumPy 中的高效 3x3 和 2x2 行列式
【发布时间】：2012-10-26 00:20:35
【问题描述】：

我有一个大的 numpy 数组 arr，形状为 (N, D, M, D)，其中 D 是两个或三个。该数组可以被认为是在N 和M 维度中被块在一起的(D,D) 矩阵块。我希望采用每个(D,D) 矩阵的行列式：

out = np.empty((N,M))
for i in xrange(N):
    for j in xrange(M):
        out[i,j] = np.linalg.det(arr[i,:,j,:])

唯一的问题是np.linalg.det 没有特殊情况的小矩阵，因此每个块的完整 BLAS 调用和 LU 分解也是如此。事实上，对于一个 2 x 2 矩阵 np.linalg.det 在我最近的 Core i7 系统上需要大约 40us。我有哪些方法可以提高该片段的性能？

【问题讨论】：

标签： python numpy linear-algebra

【解决方案1】：

为了阐明 duffymo 的观点，您可以使用向量化的numpy 操作轻松利用自己的函数利用块结构。这是否会给您带来加速可能取决于 M 和 N 的大小：

>>> a = numpy.arange(2 * 2 * 2 * 2).reshape(2, 2, 2, 2)
>>> a[:,0,:,0] * a[:,1,:,1] - a[:,1,:,0] * a[:,0,:,1]
array([[-4, -4],
       [-4, -4]])
>>> a = numpy.arange(2 * 3 * 2 * 3).reshape(2, 3, 2, 3)
>>> a[:,0,:,0] * a[:,1,:,1] * a[:,2,:,2] + \
    a[:,0,:,1] * a[:,1,:,2] * a[:,2,:,0] + \
    a[:,0,:,2] * a[:,1,:,0] * a[:,2,:,1] - \
    a[:,0,:,0] * a[:,1,:,2] * a[:,2,:,1] - \
    a[:,0,:,1] * a[:,1,:,0] * a[:,2,:,2] - \
    a[:,0,:,2] * a[:,1,:,1] * a[:,2,:,0]
array([[0, 0],
       [0, 0]])

【讨论】：

您选择 12 毫升而不是 9 毫升的任何原因？（您可以将第一项分解为 a[:,0,:,0]*(a[:,1,:,1] * a[:,2,:,2] - a[:,1,: ,2] * a[:,2,:,1]) 等）
@FreddieWitherden，没有理由——我只是懒得想了半秒多：)

【解决方案2】：

您不需要 LU 分解来计算 2x2 或 3x3 矩阵的行列式。公式是众所周知的。为什么不编写自己的函数并调用它们？

【讨论】：

两个原因：首先，我用 Python 编写的任何东西都会比 NumPy 在 C 中可能已经拥有的任何东西慢很多。其次，我也对利用块结构感兴趣（避免显式的 for循环）。
我不同意。 Python中的单行公式会比LU decomp慢吗？