加速计算矩阵辅因子的python代码答案

【问题标题】：Speed up python code for computing matrix cofactors加速计算矩阵辅因子的python代码
【发布时间】：2011-09-25 12:56:17
【问题描述】：

作为复杂任务的一部分，我需要计算 matrix cofactors。我使用nice code for computing matrix minors 以一种简单的方式做到了这一点。这是我的代码：

def matrix_cofactor(matrix):
    C = np.zeros(matrix.shape)
    nrows, ncols = C.shape
    for row in xrange(nrows):
        for col in xrange(ncols):
            minor = matrix[np.array(range(row)+range(row+1,nrows))[:,np.newaxis],
                           np.array(range(col)+range(col+1,ncols))]
            C[row, col] = (-1)**(row+col) * np.linalg.det(minor)
    return C

原来这个矩阵辅因子代码是瓶颈，我想优化一下上面的代码sn-p。关于如何做到这一点的任何想法？

【问题讨论】：

一般的bottleneck-killer是用C写瓶颈。这里有技术人员吗？
想详细说明为什么需要计算“辅因子”？是否有可能避免它并尝试为您的问题找到更直接的解决方案？即使遵循borribles 的建议，你也不会接近这样的加速，从“适当的角度”（如果可能的话）解释问题时可能会发生什么。谢谢
@eat，无法避免它们。太复杂了，这里就不解释了……
@Jeff：请注意，根据pvs 的回答，存在一个真正的障碍，你不能指望超越。因此，只有算法上的“创新”才能真正让你走得更远（至少在性能方面）。谢谢
如果您尝试计算行列式，最好的方法是关键凝聚的 chio 规则。

标签： python matrix performance numpy linear-algebra

【解决方案1】：

我建议不要使用逆行列式，而是使用 SVD

def cofactors(A):
    U,sigma,Vt = np.linalg.svd(A)
    N = len(sigma)
    g = np.tile(sigma,N)
    g[::(N+1)] = 1
    G = np.diag(-(-1)**N*np.product(np.reshape(g,(N,N)),1)) 
    return U @ G @ Vt

【讨论】：

能否请您提供此算法的参考？我对 g 和 G 发生的事情特别感兴趣。我见过基于 SVD 的算法，但它们也使用了行列式。例如：scicomp.stackexchange.com/questions/33028/…
这个推导来自对辅因子矩阵的一般定义应用单值分解。 C = det(A) A^{-T} 等于 det(U) det(S) det(V) U S^{-1} V^T。 det(U) 和 det(V) 只是 +1 或 -1，因为它们是正交矩阵，而 G 是 det(S) S^{-1} 的约简（有时标记为 Gamma）。这一切都以 det(U)det(V) U G V^{T} 的形式出现（如上所示）

【解决方案2】：

from sympy import *
A = Matrix([[1,2,0],[0,3,0],[0,7,1]])
A.adjugate().T

输出（即辅因子矩阵）为：

Matrix([
[ 3, 0,  0],
[-2, 1, -7],
[ 0, 0,  3]])

【讨论】：

【解决方案3】：

如果您的矩阵是可逆的，则辅因子与逆相关：

def matrix_cofactor(matrix):
    return np.linalg.inv(matrix).T * np.linalg.det(matrix)

这会带来很大的加速（对于 50x50 的矩阵，大约是 1000 倍）。主要原因是根本性的：这是一个O(n^3) 算法，而基于minor-det 的算法是O(n^5)。

这可能意味着对于不可逆矩阵，也有一些聪明的方法来计算辅因子（即，不使用上面使用的数学公式，而是使用其他等效定义）。

如果你坚持使用基于 det 的方法，你可以做以下事情：

大部分时间似乎都在det 内度过。（请查看 line_profiler 以自行查找。）您可以尝试通过将 Numpy 与英特尔 MKL 链接来加速该部分，但除此之外，没有什么可以做的。

你可以像这样加速其他部分的代码：

minor = np.zeros([nrows-1, ncols-1])
for row in xrange(nrows):
    for col in xrange(ncols):
        minor[:row,:col] = matrix[:row,:col]
        minor[row:,:col] = matrix[row+1:,:col]
        minor[:row,col:] = matrix[:row,col+1:]
        minor[row:,col:] = matrix[row+1:,col+1:]
        ...

这会增加大约 10-50% 的总运行时间，具体取决于矩阵的大小。原始代码有 Python range 和列表操作，它们比直接切片索引要慢。您也可以尝试更聪明，只复制实际更改的次要部分 --- 然而，在上述更改之后，接近 100% 的时间都花在 numpy.linalg.det 内，以便进一步优化其他部分没有多大意义。

【讨论】：

优秀的答案！我的矩阵是可逆的，因此一个衬里可以节省大量时间。
这会计算伴随矩阵，而不是辅因子矩阵。 det(A) * inverse(A) = adjoint(A)
@v3ga：请实际阅读答案。它计算det(A) * inverse(A)^T。辅因子是副词的转置。

【解决方案4】：

np.array(range(row)+range(row+1,nrows))[:,np.newaxis] 的计算不依赖于col，因此您可以将其移到内部循环之外并缓存该值。根据您拥有的列数，这可能会进行少量优化。

【讨论】：