Ускорение кода 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

Если ваша матрица обратима, кофактор связан с обратным:

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

Это дает большие ускорения (~ 1000x для матриц 50x50).Основная причина фундаментальна: это алгоритм 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, так что дальнейшая оптимизация других частейне имеет особого смысла.

Answer 2

Расчет np.array(range(row)+range(row+1,nrows))[:,np.newaxis] не зависит от col, так что вы можете переместить это за пределы внутреннего цикла и кэшировать значение. В зависимости от количества имеющихся у вас столбцов это может дать небольшую оптимизацию.

Answer 3

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