Inwersja macierzy bez Numpy

Question

Chcę odwrócić macierz bez użycia numpy.linalg.inv.

Powodem jest to, że używam Numba do przyspieszenia kodu, ale numpy.linalg.inv nie jest obsługiwany, więc zastanawiam się, czy mogę odwrócić macierz z "klasycznym" kodem Pythona.

Z numpy.linalg.inv przykładowy kod będzie wyglądał tak:

import numpy as np 
M = np.array([[1,0,0],[0,1,0],[0,0,1]]) 
Minv = np.linalg.inv(M) 

Answer 1

, że stosuje się wzór z http://cg.info.hiroshima-cu.ac.jp/~miyazaki/knowledge/teche23.html zapisu funkcji, które wykonuje się inwersję macierzy 4x4:

import numpy as np 

def myInverse(A): 
    detA = np.linalg.det(A) 

    b00 = A[1,1]*A[2,2]*A[3,3] + A[1,2]*A[2,3]*A[3,1] + A[1,3]*A[2,1]*A[3,2] - A[1,1]*A[2,3]*A[3,2] - A[1,2]*A[2,1]*A[3,3] - A[1,3]*A[2,2]*A[3,1] 
    b01 = A[0,1]*A[2,3]*A[3,2] + A[0,2]*A[2,1]*A[3,3] + A[0,3]*A[2,2]*A[3,1] - A[0,1]*A[2,2]*A[3,3] - A[0,2]*A[2,3]*A[3,1] - A[0,3]*A[2,1]*A[3,2] 
    b02 = A[0,1]*A[1,2]*A[3,3] + A[0,2]*A[1,3]*A[3,1] + A[0,3]*A[1,1]*A[3,2] - A[0,1]*A[1,3]*A[3,2] - A[0,2]*A[1,1]*A[3,3] - A[0,3]*A[1,2]*A[3,1] 
    b03 = A[0,1]*A[1,3]*A[2,2] + A[0,2]*A[1,1]*A[2,3] + A[0,3]*A[1,2]*A[2,1] - A[0,1]*A[1,2]*A[2,3] - A[0,2]*A[1,3]*A[2,1] - A[0,3]*A[1,1]*A[2,2] 

    b10 = A[1,0]*A[2,3]*A[3,2] + A[1,2]*A[2,0]*A[3,3] + A[1,3]*A[2,2]*A[3,0] - A[1,0]*A[2,2]*A[3,3] - A[1,2]*A[2,3]*A[3,0] - A[1,3]*A[2,0]*A[3,2] 
    b11 = A[0,0]*A[2,2]*A[3,3] + A[0,2]*A[2,3]*A[3,0] + A[0,3]*A[2,0]*A[3,2] - A[0,0]*A[2,3]*A[3,2] - A[0,2]*A[2,0]*A[3,3] - A[0,3]*A[2,2]*A[3,0] 
    b12 = A[0,0]*A[1,3]*A[3,2] + A[0,2]*A[1,0]*A[3,3] + A[0,3]*A[1,2]*A[3,0] - A[0,0]*A[1,2]*A[3,3] - A[0,2]*A[1,3]*A[3,0] - A[0,3]*A[1,0]*A[3,2] 
    b13 = A[0,0]*A[1,2]*A[2,3] + A[0,2]*A[1,3]*A[2,0] + A[0,3]*A[1,0]*A[2,2] - A[0,0]*A[1,3]*A[2,2] - A[0,2]*A[1,0]*A[2,3] - A[0,3]*A[1,2]*A[2,0] 

    b20 = A[1,0]*A[2,1]*A[3,3] + A[1,1]*A[2,3]*A[3,0] + A[1,3]*A[2,0]*A[3,1] - A[1,0]*A[2,3]*A[3,1] - A[1,1]*A[2,0]*A[3,3] - A[1,3]*A[2,1]*A[3,0] 
    b21 = A[0,0]*A[2,3]*A[3,1] + A[0,1]*A[2,0]*A[3,3] + A[0,3]*A[2,1]*A[3,0] - A[0,0]*A[2,1]*A[3,3] - A[0,1]*A[2,3]*A[3,0] - A[0,3]*A[2,0]*A[3,1] 
    b22 = A[0,0]*A[1,1]*A[3,3] + A[0,1]*A[1,3]*A[3,0] + A[0,3]*A[1,0]*A[3,1] - A[0,0]*A[1,3]*A[3,1] - A[0,1]*A[1,0]*A[3,3] - A[0,3]*A[1,1]*A[3,0] 
    b23 = A[0,0]*A[1,3]*A[2,1] + A[0,1]*A[1,0]*A[2,3] + A[0,3]*A[1,1]*A[2,0] - A[0,0]*A[1,1]*A[2,3] - A[0,1]*A[1,3]*A[2,0] - A[0,3]*A[1,0]*A[2,1] 

    b30 = A[1,0]*A[2,2]*A[3,1] + A[1,1]*A[2,0]*A[3,2] + A[1,2]*A[2,1]*A[3,0] - A[1,0]*A[2,1]*A[3,2] - A[1,1]*A[2,2]*A[3,0] - A[1,2]*A[2,0]*A[3,1] 
    b31 = A[0,0]*A[2,1]*A[3,2] + A[0,1]*A[2,2]*A[3,0] + A[0,2]*A[2,0]*A[3,1] - A[0,0]*A[2,2]*A[3,1] - A[0,1]*A[2,0]*A[3,2] - A[0,2]*A[2,1]*A[3,0] 
    b32 = A[0,0]*A[1,2]*A[3,1] + A[0,1]*A[1,0]*A[3,2] + A[0,2]*A[1,1]*A[3,0] - A[0,0]*A[1,1]*A[3,2] - A[0,1]*A[1,2]*A[3,0] - A[0,2]*A[1,0]*A[3,1] 
    b33 = A[0,0]*A[1,1]*A[2,2] + A[0,1]*A[1,2]*A[2,0] + A[0,2]*A[1,0]*A[2,1] - A[0,0]*A[1,2]*A[2,1] - A[0,1]*A[1,0]*A[2,2] - A[0,2]*A[1,1]*A[2,0] 

    Ainv = np.array([[b00, b01, b02, b03], [b10, b11, b12, b13], [b20, b21, b22, b23], [b30, b31, b32, b33]])/detA 

return Ainv 

Answer 2

Przez 4 x 4 matrycy to chyba tylko o OK, aby użyć formuły matematyczne, które można znaleźć za pomocą googlowanie "formuła dla odwrotności macierzy 4 na 4". Na przykład tutaj (nie może zagwarantować dokładność)

http://www.cg.info.hiroshima-cu.ac.jp/~miyazaki/knowledge/teche23.html

Ogólnie wywracania zwykłą macierz nie jest słaby serca. Musisz być świadomy wszystkich matematycznie trudnych przypadków i wiedzieć, dlaczego nie będą one miały zastosowania do twojego użycia i złapać je, gdy otrzymasz matematycznie patologiczne dane wejściowe (które, lub zwróć wyniki o niskiej dokładności lub numeryczne śmieci w wiedzy, że nie będzie miało to znaczenia w przypadku użycia, pod warunkiem, że w rzeczywistości nie dojdzie do podziału przez zero lub przepełnienia MAXFLOAT ..., który można złapać za pomocą obsługi wyjątku i przedstawić jako "Błąd: macierz jest pojedyncza lub bardzo blisko".

Generalnie lepiej jako programista użyć kodu bibliotecznego napisanego przez ekspertów matematyki numerycznej, chyba że jesteś gotów poświęcić czas na zrozumienie fizycznej i matematycznej natury konkretnego problemu, który rozwiązujesz i stajesz się ekspertem w dziedzinie matematyki we własnym zakresie specjalistyczna dziedzina.

Answer 3

tutaj jest bardziej elegancka i skalowalne rozwiązanie MOM. Będzie działać dla każdej macierzy nxn i możesz znaleźć zastosowanie dla innych metod. Zauważ, że getMatrixInverse (m) przyjmuje tablicę tablic jako dane wejściowe. Zachęcamy do zadawania pytań.

def transposeMatrix(m): 
    return map(list,zip(*m)) 

def getMatrixMinor(m,i,j): 
    return [row[:j] + row[j+1:] for row in (m[:i]+m[i+1:])] 

def getMatrixDeternminant(m): 
    #base case for 2x2 matrix 
    if len(m) == 2: 
     return m[0][0]*m[1][1]-m[0][1]*m[1][0] 

    determinant = 0 
    for c in range(len(m)): 
     determinant += ((-1)**c)*m[0][c]*getMatrixDeternminant(getMatrixMinor(m,0,c)) 
    return determinant 

def getMatrixInverse(m): 
    determinant = getMatrixDeternminant(m) 
    #special case for 2x2 matrix: 
    if len(m) == 2: 
     return [[m[1][1]/determinant, -1*m[0][1]/determinant], 
       [-1*m[1][0]/determinant, m[0][0]/determinant]] 

    #find matrix of cofactors 
    cofactors = [] 
    for r in range(len(m)): 
     cofactorRow = [] 
     for c in range(len(m)): 
      minor = getMatrixMinor(m,r,c) 
      cofactorRow.append(((-1)**(r+c)) * getMatrixDeternminant(minor)) 
     cofactors.append(cofactorRow) 
    cofactors = transposeMatrix(cofactors) 
    for r in range(len(cofactors)): 
     for c in range(len(cofactors)): 
      cofactors[r][c] = cofactors[r][c]/determinant 
    return cofactors