Przyspieszenie kodu Pythona do obliczania kofaktorów macierzy

Question

W ramach złożonego zadania muszę obliczyć matrix cofactors. Zrobiłem to w prosty sposób, używając tego nice code for computing matrix minors. Tu jest mój kodu:

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 

Okazuje się, że ten kod kofaktorem matryca jest wąskim gardłem, i chciałbym, aby zoptymalizować fragment kodu powyżej. Wszelkie pomysły, jak to zrobić?

Answer 1

Jeśli macierz jest odwracalna, kofaktor związany jest z odwrotności:

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

To daje duże speedups (~ 1000x do 50x50 matryce). Główny powód jest fundamentalny: jest to algorytm O(n^3), podczas gdy ten oparty na pomniejszym detekcie to O(n^5).

To prawdopodobnie oznacza, że ​​również dla nieuszkodzonych matryc istnieje pewien sprytny sposób obliczenia kofaktora (tj. Nie używaj wzoru matematycznego, którego używasz powyżej, ale jakiejś innej równoważnej definicji).

---

Jeśli trzymać się z podejścia opartego det, co można zrobić, jest następujący:

Większość czasu wydaje się być wydane wewnątrz det. (Aby to sprawdzić, sprawdź: line_profiler) Możesz spróbować przyspieszyć tę część, łącząc Numpy z Intel MKL, ale poza tym niewiele można zrobić.

Można przyspieszyć drugą część kodu tak:

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

ten zyskuje jakieś 10-50% całkowitego czasu pracy w zależności od wielkości firmy matryc. Oryginalny kod ma Python range i manipuluje listami, które są wolniejsze niż bezpośrednie indeksowanie plasterków. Możesz też spróbować być bardziej sprytnym i skopiować tylko części mniejszego, które faktycznie się zmieniają - jednak już po powyższej zmianie prawie 100% czasu spędzasz wewnątrz numpy.linalg.det, więc dodatkowa optymalizacja pozostałych części nie mieć taki sens.

Answer 2

Obliczenie np.array(range(row)+range(row+1,nrows))[:,np.newaxis] nie zależało od col, więc mógłbyś przenieść to poza wewnętrzną pętlę i zapisać wartość w pamięci podręcznej. W zależności od liczby kolumn może to dać małą optymalizację.