Inversion de matrice sans Numpy

Je veux inverser une matrice sans l'aide de numpy.linalg.inv.

La raison en est que je suis en utilisant Numba pour accélérer le code, mais numpy.linalg.inv n'est pas pris en charge, donc je me demande si je peux inverser une matrice avec des "classiques" du code Python.

Avec numpy.linalg.inv un exemple de code devrait ressembler à ça:

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

Probablement pas. Il n'y a pas de python "builtin" le faire pour vous et la programmation d'une inversion de matrice vous-même est tout sauf facile (voir par exemple en.wikipedia.org/wiki/... pour une liste (probablement non exhaustive de méthodes). Je suis pas au courant de tout numpyindépendant de package d'algèbre linéaire pour python...
Si vous voulez inverser des matrices 3x3 seulement, vous pouvez consulter la formule ici. (Il vaut mieux spécifier la dimension et le type de matrices que vous souhaitez inverser. Dans votre exemple vous utilisez le plus trivial matrice d'identité. Sont-ils réels? Et régulier?)
Pour être précis, est une véritable matrice 4x4

OriginalL'auteur Alessandro Vianello | 2015-08-20

31

Ici est un plus élégant et d'une solution évolutive, de l'omi. Il va travailler pour toute matrice nxn et vous peuvent trouver une utilisation pour les autres méthodes. Notez que getMatrixInverse(m) dans un tableau de tableaux en entrée. N'hésitez pas à poser des questions.
```
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
```
Cela fonctionne parfaitement. En fonction de l'exigence, doit être la accepté de répondre. Le seul petit changement nécessaire est en #base case for 2x2 matrix. vous avez besoin de convertir explicitement à flotteur.
Si la matrice n'est pas carrée la fonction de transposition donnera une erreur, pour trouver le transposer pour une liste simplement que nous pouvons faire: zip(*theArray) Prises à partir de: stackoverflow.com/questions/4937491/matrix-transpose-in-python
vous avez raison, merci. Bien que non matrices carrées n'ont pas l'inverse, je ne demande ma réponse est composée de pièces réutilisables, donc j'ai corrigé la fonction de transposition à votre suggestion.
J'ai utilisé une partie de votre code avec une matrice non carrée pour trouver la matrice de projection, c'est pourquoi j'ai eu cette question. Merci
c'est formidable. Je souhaite que je pourrais upvote plus d'une fois

OriginalL'auteur stackPusher
3

Pour un 4 x 4 de la matrice, il est probablement juste sur OK pour utiliser la formule mathématique, que vous pouvez trouver à l'aide de Googler "formule de 4 par 4 de la matrice inverse". Par exemple ici (je ne peux pas se porter garant de son exactitude):

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

En général, l'inversion d'une matrice n'est pas pour les timides. Vous devez être conscient de tous les mathématiquement les cas difficiles et de savoir pourquoi ils ne s'appliquent à votre utilisation, et les attraper quand vous êtes fourni avec mathématiquement pathologique entrées (ou revenir aux résultats de la faible précision ou numérique ordures en sachant qu'il ne sera pas question dans votre cas d'utilisation à condition de ne pas parvenir à diviser par zéro ou de débordement MAXFLOAT ... que vous pourriez attraper avec un gestionnaire d'exception et de les présenter comme "Erreur: la matrice est singulière ou très proche de celle-ci").

Il est généralement meilleur en tant que programmeur pour utiliser la bibliothèque de code écrit par le numérique en mathématiques experts, sauf si vous êtes prêt à passer du temps à la compréhension de la physique et de la mathématique, de la nature du problème particulier qui vous vous adressez et de devenir votre propre mathématiques expert dans votre domaine d'expertise.

OriginalL'auteur nigel222

À tout le moins au 16 juillet 2018 Numba a un rapide de la matrice inverse. (Vous pouvez voir comment ils surcharge de la norme NumPy inverse et d'autres opérations ici.)

Voici les résultats de mon analyse comparative:

import numpy as np
from scipy import linalg as sla
from scipy import linalg as nla
import numba
def gen_ex(d0):
x = np.random.randn(d0,d0)
return x.T + x
@numba.jit
def inv_nla_jit(A):
return np.linalg.inv(A)
@numba.jit
def inv_sla_jit(A):
return sla.inv(A)

Pour les petites matrices, il est particulièrement rapide:

ex1 = gen_ex(4)
%timeit inv_nla_jit(ex1) # NumPy + Numba
%timeit inv_sla_jit(ex1) # SciPy + Numba
%timeit nla.inv(ex1)     # NumPy
%timeit sla.inv(ex1)     # SciPy

[Out]

2.54 µs ± 467 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
67.3 µs ± 9.18 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
63.5 µs ± 7.65 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
56.6 µs ± 5.03 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

Avis que l'accélération ne fonctionne que pour NumPy inverse, pas SciPy (comme prévu).

Légèrement plus grande de la matrice:

ex2 = gen_ex(40)
%timeit inv_nla_jit(ex2) # NumPy + Numba
%timeit inv_sla_jit(ex2) # SciPy + Numba
%timeit nla.inv(ex2)     # NumPy
%timeit sla.inv(ex2)     # SciPy

[Out]

131 µs ± 12.9 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
278 µs ± 26.2 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
231 µs ± 24.5 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
189 µs ± 11.2 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

Donc, il y a toujours un speedup ici, mais SciPy est en train de rattraper.

OriginalL'auteur webelo

-5

J'ai utilisé la formule de http://cg.info.hiroshima-cu.ac.jp/~miyazaki/knowledge/teche23.html d'écrire la fonction qui fait l'inversion d'une 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

Vous ne voulez pas utiliser np.linalg.inv mais np.linalg.det est bien? C'est une très maladroit exigence...
Bien sûr, on a besoin d'écrire un autre "brute force" mise en œuvre pour le déterminant de calcul en tant que bien. Ou simplement calculer le det à l'extérieur de la Numba de la fonction et de la passer en argument

OriginalL'auteur Alessandro Vianello

Vous devez vous connecter pour publier un commentaire.