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

Question

在一个复杂的任务中，我需要计算矩阵的余子式。我用一种简单的方法实现了这一点，使用了这个计算矩阵小行列式的好代码。以下是我的代码：

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

结果发现，这段计算矩阵余子式的代码是性能瓶颈，我想对上面的代码进行优化。有没有什么好的建议呢？

Answer 1

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

Answer 2

这个计算 np.array(range(row)+range(row+1,nrows))[:,np.newaxis] 的结果和 col 没有关系，所以你可以把这个计算放到内层循环外面，这样可以把结果存起来。根据你有多少列，这样做可能会稍微提高一点效率。

Answer 3

如果你的矩阵是可逆的，那么余子式和逆矩阵是有关系的：

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

这样做可以大大加快速度（对于50x50的矩阵，速度提升大约是1000倍）。主要原因是这个算法的基本性质：这是一个O(n^3)的算法，而基于小行列式的算法是O(n^5)。

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

---

如果你继续使用基于行列式的方法，你可以这样做：

大部分时间似乎都花在了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里，所以进一步优化其他部分就没什么意义了。