Acelera el código python para calcular cofactores de matrices

Question

Como parte de una tarea compleja, necesito calcular matrix cofactors. Hice esto de manera directa usando nice code for computing matrix minors. Aquí está mi código:

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 

Resulta que este código matriz cofactor es el cuello de botella, y me gustaría para optimizar el fragmento de código anterior. ¿Alguna idea de cómo hacer esto?

Answer 1

Si su matriz es invertible, el cofactor está relacionado con el inverso:

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

Esto da grandes aceleraciones (~ 1000x de 50x50 matrices). La razón principal es fundamental: este es un algoritmo O(n^3), mientras que el menor basado en det es O(n^5).

Esto probablemente significa que también para las matrices no invertibles, hay alguna forma inteligente de calcular el cofactor (es decir, no usar la fórmula matemática que usted usa arriba, sino alguna otra definición equivalente).

---

Si usted se pega con el enfoque basado en det, lo que puede hacer es la siguiente:

La mayoría de las veces parece ser gastado en el interior det. (Consulte line_profiler para averiguarlo usted mismo.) Puede intentar acelerar esa parte vinculando Numpy con Intel MKL, pero aparte de eso, no hay mucho que se pueda hacer.

Puede acelerar la otra parte del código como este:

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:] 
     ... 

Este gana un poco de tiempo de ejecución total del 10-50%, dependiendo del tamaño de sus matrices. El código original tiene Python range y listas de manipulaciones, que son más lentas que la indexación de sectores directos. También podría tratar de ser más astuto y copiar solo las partes del menor que realmente cambian, sin embargo, ya después del cambio anterior, cerca del 100% del tiempo se usa dentro de numpy.linalg.det para que la optimización de las otras partes no lo haga tiene mucho sentido.

Answer 2

El cálculo de np.array(range(row)+range(row+1,nrows))[:,np.newaxis] no depende de col por lo que podría moverlo fuera del bucle interno y guardar el valor en caché. Dependiendo de la cantidad de columnas que tenga, esto podría dar una pequeña optimización.