Inversion numérique de (très) grandes matrices

[Anakinele]

Bonjour,

Dans un programme en Python (mais peu importe le langage en fait), je suis amené à inverser une matrice de très grande taille (facilement plusieurs millions de cases, beaucoup plus si la RAM le permettait) contenant des coefficients variant au plus sur 2-3 ordres de grandeurs et positifs. Seuls les coefficients diagonaux sont toujours exactement nuls et, enfin, la matrice est symétrique.

Comme je pouvais m'y attendre, l'inversion de telles matrices me pose des problèmes...
Le premier d'entres-eux (peut-être le seul mais je n'ai pas réussi à aller plus loin pour l'instant) est que le déterminant calculé numériquement est soit nul soit infini. Ca ne m'étonne pas mais ça m'ennuie beaucoup...

Je me suis dit qu'en normalisant la matrice avec quelque chose du genre det(M/max(M)) je pourrais m'en sortir mais ça ne marche que dans certains cas... J'ai aussi essayé de normaliser, puis de calculer le déterminant puis -- avec une inspiration soudaine -- de "normaliser" une seconde fois en divisant par la racine n-ème du déterminant de la matrice normalisée avec n l'ordre de ma matrice. Bref, j'ai tatonné un peu n'importe comment, j'ai trouvé un tas de solutions bancales pour des cas particuliers mais en 30 secondes j'arrive toujours à ressortir une nouvelle matrice dont je ne peux pas calculer correctement le déterminant numériquement... la finalité étant, pour rappel, d'inverser la matrice. Sinon, je pourrais utiliser un algorithme d'inversion ne nécessitant pas de calculer le déterminant, je ne sais pas comment fonction inv() de Python, je vais allez voir ça...

Merci d'avance pour toute aide ou piste de réflexion!
-----

[toothpick-charlie]

le logiciel R est assez bon pour inverser des matrices. Comme c'est de l'open-source, tu peux obtenir le code. Ou sinon tu peux appeler R depuis ton code.

[Tryss]

L'inversion de matrices est un des principaux problèmes de calcul numérique.

Truc important : ne jamais calculer le déterminant, c'est aussi couteux que l'inversion en elle même.


Après il y a tout un tas de méthodes différentes, qui dépendent de ta matrice (par exemple est elle définie positive?)

[obi76]

Bonjour,

si vos matrices sont vraiment très très grandes, pour les notres (50000x50000) qui sont relativement creuses, on utilise une bibliothèque développée par l'INRIA nommée PASTIX.

[A voir en vidéo sur Futura]

[Anakinele]

Merci à tous pour vos réponses.

@ toothpick-charlie et obi76 :
Je ne suis pas contre utiliser un programme externe à Python mais si je pars dans cette direction, je pense que j'aurai l'embarras du choix. Il y a Lapack aussi... Je me rends compte que j'avais sous-estimé la difficulté numérique d'inversion de très grandes matrices mais a posteriori ça m'aurait arrangé que Python gère ça en interne... Pour d'autres fonctions, je sais que Python s'appuie sur des codes de bibliothèques robustes en Fortran ou C, si c'est le cas pour la fonction inv() d'inversion de matrices, je ne suis pas sûr que l'utilisation d'un autre langage/programme résolve le problème facilement, mais je vais me pencher sur la question. Le problème est surtout -- comme toujours -- de faire communiquer les différents codes... j'appréhende déjà

En ce qui concerne PASTIX, je garde ça en tête, sait-on jamais, mais ici ma matrice est loin d'être creuse (seuls les termes diagonaux sont nuls), sinon pour les matrices creuses, je sais qu'il existe de nombreux algorithme astucieux facilitant les calculs ; Python en propose visiblement.

@ Tryss :
Je sais que ce n'était qu'une remarque mais le calcul du déterminant me permettait juste de vérifier que Python ait une chance de pouvoir inverser la matrice. Peut-être que l'algorithme Python utilisé ne se préoccupe pas a priori de la valeur du déterminant mais pour l'instant j'ai toujours observé que quand le déterminant est numériquement nul ou infini, l'inversion échoue.
Pour la nature exacte de mes matrices, je ne peux pas être plus précis que ce que j'ai indiqué dans mon premier post. Elles sont calculées numériquement et présentent du coup peu de propriétés évidentes.

[invitefb0f1e11]

Salut,

Je te conseille de faire une SVD, comme ça tu as accès aux valeurs singulière et donc tu sait directement si ta matrice est singulière où non (présence de valeurs singulière nulle).
Ensuite tu n'as plus qu'à inversé/pseudo-inversé (si la matrice est singulière) en prenant l'inverse de de la matrice diagonale et en recomposant la matrice. En oubliant pas que :



Les matrice U et V étant unitaire.

Je viens de faire des test avec la librairie LAPACK, pour des matrice 5000x5000 (plusieurs millions d'éléments) et le temps de calcule est "raisonnable".
Sinon il me semble que SciPy intègre une partie des routine LAPACK ... à voir.

@+,
G.

[toothpick-charlie]

en tout cas à ta place j'installerais R sur ma bécane (c'est gratos) et j'essaierai d'inverser quelques-unes des matrices qui te donnent du fil à retordre (c'est la fonction "solve"). Chez moi une matrice 1000x1000 est inversée en 1 seconde environ (I5 2.66GHz)

[Anakinele]

@ toothpick-charlie :
Bon, j'essaierai R alors, depuis le temps que j'en ai entendu parler. La curiosité est déjà une bonne raison
Ceci dit, une matrice 1000x1000 ne me pose pas toujours de problème (ça dépend des coefficients de la matrice) et s'inverse aussi très vite en Python. Ca commence à devenir plus tendu, quels que soient les coefficients pour une matrice 5000x5000 et à 10000x10000 c'est plus possible à cause de la RAM qui dit "ouf"

@ Psyricien :
J'ai essayé SVD -- c'est que je n'y avais pas pensé -- mais en fait, j'ai aussi un problème de taille de matrice pour la RAM et cette solution ne résolverait pas tout.


Finalement, j'ai écrit un code récursif d'inversion de matrices par bloc pour inverser des matrices de taille arbitraire écrite sur le disque dur (avec h5py pour ceux qui connaissent). J'ai ainsi pu inverser une matrice 10000x10000 (sans coefficient nuls, soit 800Mo de données) sans utiliser plus de quelques dizaines de Mo de RAM en environ 20 minutes sur un bon PC (sans plus !). Numpy ne permettait pas de faire ça avec inv()... pas en 32 bits avec mes 2-3Go de RAM disponibles ! Le code tourne bien, j'ai juste un problème de précision qui s'accroit avec la taille des matrices (mais ça devrait me suffire).

Je testerai d'autres solutions mais à court terme je me contenterai de ça pour ne pas me prendre la tête à interfacer Python avec un autre langage !

Merci de vos participations !
Bon WE