Inversione matrice senza Numpy

Question

Voglio invertire una matrice senza utilizzare numpy.linalg.inv.

Il motivo è che sto usando Numba per accelerare il codice, ma numpy.linalg.inv non è supportato, quindi mi chiedo se posso invertire una matrice con codice Python "classico".

Con numpy.linalg.inv un codice di esempio che sarebbe simile:

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

Answer 1

ho usato la formula dalla http://cg.info.hiroshima-cu.ac.jp/~miyazaki/knowledge/teche23.html per scrivere la funzione che fa l'inversione di una matrice 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

Per una matrice 4 x 4 è probabilmente solo su OK per utilizzare la formula matematica, che potete trovare tramite Googling "formula per 4 per 4 matrice inversa". Per esempio qui (non posso garantire per la precisione):

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

In generale invertendo una matrice generale non è per i deboli di cuore. Devi essere consapevole di tutti i casi matematicamente difficili e sapere perché non si applicheranno al tuo utilizzo, e prenderli quando ti verranno forniti input matematici patologici (o restituire risultati di bassa precisione o immondizia numerica nella consapevolezza che non importa nel tuo caso d'uso a patto che tu non finisca per dividere per zero o traboccare MAXFLOAT ... che potresti prendere con un gestore di eccezioni e presentare come "Errore: matrice è singolare o molto vicina ad esso").

In generale, è meglio che un programmatore utilizzi il codice di libreria scritto da esperti di matematica numerica, a meno che non si sia disposti a dedicare del tempo a comprendere la natura fisica e matematica del problema specifico che si sta affrontando e diventare il proprio esperto di matematica nel proprio campo specialistico.

Answer 3

Ecco una soluzione più elegante e scalabile, imo. Funzionerà per qualsiasi matrice nxn e potresti trovare uso per gli altri metodi. Nota che getMatrixInverse (m) prende in input una matrice di matrici. Non esitate a fare domande.

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