加速Python代码以计算矩阵余子式

Question

加速Python代码以计算矩阵余子式

pythonmatrixperformancenumpylinear-algebra

10

作为一个复杂任务的一部分，我需要计算矩阵余子式。我使用了计算矩阵子式的好代码来简单地完成这个任务。以下是我的代码：

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

事实证明，矩阵余子式代码是瓶颈所在，我想优化上面的代码片段。有什么想法吗？

- new name

通用的瓶颈克星是使用 C 语言编写瓶颈代码。这里有技术人员吗？ - Evpok

能否详细说明为什么需要计算“余子式”？是否可能避免它并尝试找到更直接的解决方案来解决问题？即使按照“borrible”的建议，您也不会获得任何接近这样的加速，除非从“正确的角度”（如果可能）解释问题。谢谢。 - eat

@eat，无法避免它们。在这里解释起来太复杂了... - new name

请注意，根据 pv 的回答，存在一个你不能期望越过的真实障碍。因此，只有算法上的“创新”才是真正使你能够在性能方面实现重大进展的手段。谢谢。 - eat

如果您尝试计算行列式，最好的方法是使用基本极化消元法。 - Robert William Hanks

在数值线性代数中计算行列式的常规方法是计算LU分解并从三角形L和U中找到行列式。det和LU的计算复杂度基本相同，由于已经存在针对LU的优化库，因此很难超越它们。 - pv.

4个回答

2

np.array(range（row）+range（row+1，nrows）)[:,np.newaxis]的计算与col无关，因此您可以将其移到内循环之外并缓存该值。根据您拥有的列数，这可能会带来小的优化。

- borrible

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]])

- Prokhozhii

2

建议使用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

- Greg von Winckel

请问您能提供一下这个算法的参考资料吗？我对g和G发生了什么特别感兴趣。我看过基于SVD的算法，但它们也使用了行列式。例如：https://scicomp.stackexchange.com/questions/33028/fast-algorithm-for-computing-cofactor-matrix - user2393987

这个导出式源于将单值分解应用于余子式矩阵的一般定义。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)UGV^{T}（如上所示）。 - AlphaBetaGamma96

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- pv. · Accepted Answer

如果您的矩阵可逆，那么余子式与其逆相关：

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

这大大提高了速度（对于50x50的矩阵，可达到1000倍左右）。主要原因是基本的：这是一个 O(n^3) 的算法，而基于小行列式的算法则是 O(n^5) 的。

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

如果您坚持使用基于det的方法，您可以执行以下操作：

大部分时间似乎都花费在det内部。（使用line_profiler自己找出来。）您可以尝试通过将Numpy与Intel 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内部，因此对其他部分进行进一步优化并没有太大意义。