2015-08-20 22 views
9

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

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) 
+1

Probabilmente no. Non c'è python "builtin" che lo faccia per te e programmare l'inversione di una matrice te stesso è tutt'altro che facile (vedi per esempio https://en.wikipedia.org/wiki/Invertible_matrix#Methods_of_matrix_inversion per un elenco di metodi probabilmente incompleto). Non sono a conoscenza di alcun pacchetto di algebra lineare 'numpy' indipendente per Python ... – sebastian

+0

Se vuoi invertire solo le matrici 3x3, puoi cercare la formula [qui] (http://mathworld.wolfram.com /MatrixInverse.html). (È meglio specificare la dimensione e il tipo di matrici che si desidera invertire.Nell'esempio si utilizza la matrice di identità più banale.Sono reali? E regolari?) – Falko

+0

Per essere precisi è una matrice reale 4x4 –

risposta

-4

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 
+7

Non si desidera usa 'np.linalg.inv' ma' np.linalg.det' va bene? Questo è un requisito davvero imbarazzante ... – sebastian

+0

Naturalmente è necessario scrivere un'altra implementazione "forza bruta" anche per il calcolo dei determinanti. Oppure basta calcolare il det al di fuori della funzione Numba e passarlo come argomento –

1

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.

22

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 
+0

Questo funziona perfettamente. Secondo il requisito, dovrebbe essere la risposta accettata. L'unica piccola modifica richiesta è in "#base case for 2x2 matrix'. devi convertire esplicitamente in float. – CoderFrom94

+1

Se la matrice non è quadrata, la funzione trasposizione darà un errore, per trovare la trasposizione per un elenco semplicemente possiamo fare: zip (* theArray) Tratto da: https: // stackoverflow.it/questions/4937491/matrix-transpose-in-python –

+1

@MohanadKaleia hai ragione, grazie. Sebbene le matrici non quadrate non abbiano inversioni, ritengo che la mia risposta sia composta da pezzi riutilizzabili, quindi ho corretto la funzione di trasposizione come da suggerimento. – stackPusher