Je développe actuellement un programme (en C) dans lequel j'ai besoin d'inverser une matrice. Je parle bien d'un calcul d'inverse et non de la résolution d'un système linéaire. En effet il est plus intéressant pour mon programme d'inverser une fois ma matrice pour ensuite faire des produits matriciels.
Ma matrice est carrée, peut être de taille très grande (quelques milliers **2) et n'est pas creuse. Elle est juste symétrique. C'est la matrice de covariance issue des méthodes de krigeage.
Dans la première version du programme, écrite avec R cette fois, par simplicité je résolvais le système linéaire à chaque fois, et j'ai été confronté au problème de conditonnement quand le taille augmente. J'avais réglé le problème avec un préconditionneur.
Je voudrais à présent coder l'inversion moi même en C. J'ai trouvé beaucoup de choses sur la résolution de systèmes, mais très peu sur l'inversion.
Avez-vous des infos/conseils sur une méthode qui pourrait convenir ? (Si cette méthode est parallélisable, c'est magnifique.)
Bonjour ! J'ai peur qu'il y ait peu de personnes pouvant t'aider tant le sujet est pointu...
Si le but est d'inverser la matrice et non de résoudre un système linéaire, la méthode du pivot est faite pour ça. Mais comme elle est archi connue, j'imagine que tu l'as déjà exclue. Trop lente ? Pas parallélisable ? Pour trouver plus rapide, donc moins précis, il y a les méthodes itératives, qui ont des chances d'être parallélisables, mais c'est en effet pour résoudre des systèmes linéaires. Ceci dit, inverser une matrice, c'est presque pareil, non ?
(Au fait, la matrice est symétrique, OK, mais est-elle définie positive ? ou à diagonale dominante ? Je sais qu'on peut avoir des méthodes plus efficaces dans ce genre de cas.)
Ta matrice étant symétrique, tu peux éventuellement aller jeter un oeil du côté de la décomposition LU.
C'est une décomposition qui consiste à écrire ta matrice de départ comme le produit d'une matrice triangulaire inférieure (L pour lower) et d'une matrice triangulaire supérieure (U pour upper).
J'en mettrai pas ma main à couper, mais calculer l'inverse de deux matrices avec presque la moitié des coefficients nuls et bien placés, ça devrait être plus rapide et plus facile que l'inverse d'une seule matrice complète.
Je te lance juste une piste ici, je ne sais pas si ça t'aidera ou pas...
Merci pour ta réponse, je ne manquerai pas d'y jeter un coup d'oeil. J'ai également vu des choses sur l'algorithme de Gauss Jordan qui pourrait convenir. J'essaierai de faire un petit retour ici si l'une de ces méthodes fonctionnent.
Pour les méthodes directes, avec une matrice symétrique définie positive,on utilise plutôt la factorisation de Cholesky, plus performante que LU dans ce cas . Il n'y a plus qu'une seule matrice triangulaire à déterminer avec \(A=LL^T\)
Après, autre remarque, LU c'est du \(O(n^3 )\), Cholesky du \(O(n^3/3 )\) je crois ( on gagne un peu), mais pour des très grandes tailles avec des matrices pleines , le temps de calcul peut devenir prohibitif et les méthodes itératives peut-être plus adaptées. ( mais Perso je n'ai qu'une expérience qui date déjà de matrices bandes pour des calculs de mécanique des structures par éléments finis)
- Edité par Sennacherib 28 février 2018 à 14:23:00
tout ce qui est simple est faux, tout ce qui est compliqué est inutilisable
Si j'ai bien compris, il s'agit de trouver l'inverse de la matrice car on va résoudre des tas de systèmes linéaires ayant la même matrice (c'est bien ça ?). Donc plutôt que de faire des tas de résolutions, on va juste multiplier la matrice inverse par les tas de second membres.
Si c'est bien ça le contexte, il me semble que c'est une bonne idée de réaliser une décomposition LU (ou LLt si jamais la matrice était définie positive, mais Xamime a oublié de répondre à ma question), c'est précisément fait pour cette situation. Certes, ça demande O(n³) calculs, mais juste une seule fois. Ensuite, quand on va résoudre le tas de systèmes linéaires, on n'aura presque plus rien à faire. Par contre, ce « juste une fois » est coûteux, surtout s'il n'est pas parallélisable (ce que je crains). Les méthodes itératives sont tentantes (et il me semble qu'on doit pouvoir les rendre parallélisables), mais il me semble qu'elles ne permettent pas de capitaliser sur le fait que c'est toujours le même système linéaire qu'on doit résoudre.
@robun : en effet j'ai oublié d'y répondre :/ à vrai dire, en cliquant sur la notification je n'ai vu qu'un message, je l'ai donc complétement zappé sur mon petit écran désolé
Ma matrice n'a pas d'autre propriété, ni symétrique définie positive, ni à diagonale (semi)dominante.
La factorisation LU est en effet déjà beaucoup utilisée, mais je compte tout de même l'étudier. Je ne connais pas bien les méthodes itératives et il me semblait que leur potentiel de parallélisation était assez faible, et il faudrait résoudre n (taille de la matrice) systèmes linéaires pour inverser la matrice : je les avais donc exclues. Je vais les reconsidérer. Après tout, si elles sont efficaces, je peux résoudre mes n systèmes en parallèle...
Inverser une matrice et résoude un système linéaire peut être très similaire, mais il me semble qu'il y a plus efficace. N'étant (pas encore ) expert dans le domaine, je tatonne.
"Si j'ai bien compris, il s'agit de trouver l'inverse de la matrice car on va résoudre des tas de systèmes linéaires ayant la même matrice"
C'est tout à fait ça. Dans la littérature on trouve plusieurs chiffres, et l'inversion de cette matrice prendrait jsuqu'à 80% du temps d'exécution de l'algorithme.
Pour donner davantage de détails, la matrice est construite à partir de la distance entre les points de mon échantillons de données. Les propriétés de la matrice (et notamment son conditionnement) vont dépendre de la répartition des points dans mon domaine d'étude. D'ailleurs il se peut, si mon échantillon est mal construit, que ma matrice se retrouve avec une ligne/colonne nulle, ou deux lignes/colonnes quasi identiques. D'où mon problème de conditionnement.
s'il en est encore temps, je vais apporter ma petite contribution (je pratique pas mal sur cette thématique):
En premier lieu, si tu fais du C et que tu veux comparer tes résultats avec une implémentation de référence, tu voudras utiliser les bibliothèques BLAS et LAPACK avec son wrapper en C (ce sont LES bibliothèques de référence en algèbre linéaire numérique, je pense que c'est ce qui est utilisé par R et autres). Des versions parallèles existent (OpenMP, MPI, et GPU), mais il y a "paralléliser", et "paralléliser intelligemment" ...
La seule condition pour faire un LU est que ta matrice soit inversible, elle peut donc être complexe non-symétrique etc. LAPACK te permet de faire ça (toujours si tu veux comparer).
Les méthodes itératives se parallélisent généralement très bien, tout dépend de ce que tu entends par "paralléliser" ! Elles sont basées sur des produits matrice-vecteur qui se parallélisent très bien, mais tu peux aussi lancer une inversion sur chaque thread d'exécution. Leur gros défaut est bien comme tu as pu le constater qu'on ne peut résoudre qu'un seul problème à la fois. Certaines méthodes récentes (j'ai vu passer un papier récemment...) permettent de traiter plusieurs seconds membres. Quoiqu'il en soit, si ton conditionnement est pourri, tu ne t'en sortiras pas sans pré-conditionneur. Les méthodes les plus connues sont le Gradient Conjugué ou encore l'algorithme GMRES (beaucoup plus général).
De mon expérience, la décomposition LU d'une matrice de taille N x N avec N~quelques milliers est faisable sur un seul thread tant que le "quelques milliers" ne dépasse pas ~10000 (en fait, tout dépend du temps que tu as devant toi !). Pour N=1000, ça te prendra quelques secondes, et encore...
Si tu as une ligne ou une colonne nulle, ta matrice n'est PAS inversible, et tu voudras calculer sa "pseudo-inverse" ou ajouter une condition supplémentaire.
L'inversion de la matrice est généralement ce qui prend le plus de temps dans la plupart des codes de calcul que j'ai pu voir... C'est un vrai problème sur lequel travaillent de nombreux chercheurs. Si ta matrice est par exemple construite à partir d'un potentiel gravitationnel etc. (peut-être même ta matrice de krigeage !), les méthodes les plus récentes utilisent des représentations par matrices hiérarchiques (H-matrices). Ça fonctionne très bien mais ça demande un peu d'expérience pour être implémenté.
Tu peux enfin regarder s'il n'existe pas des méthodes spécifiques à tes problèmes.
Conclusion : fais du LU ou LLt, et n'hésite pas si tu as d'autres interrogations.
Bonne journée
Avez-vous entendu parler de Julia ? Laissez-vous tenter ...
s'il en est encore temps, je vais apporter ma petite contribution (je pratique pas mal sur cette thématique):
Oui il est encore temps, je ne suis pas prêt d'avoir fini sur le sujet Merci pour ton message.
Le but est en effet de paralléliser "intelligemment". Ce n'est pas un exercice de pratique, c'est vraiment une recherche de performance. Pour le moment j'utilise GSL, mais je compte passer à ScaLapack par la suite.
Nozio a écrit:
Certaines méthodes récentes (j'ai vu passer un papier récemment...) permettent de traiter plusieurs seconds membres.
les méthodes les plus récentes utilisent des représentations par matrices hiérarchiques (H-matrices). Ça fonctionne très bien mais ça demande un peu d'expérience pour être implémenté
Si tu as la référence du papier, je suis preneur Je vais regarder du côté matrices hiérarchiques.
La liste des choses à étudier s'allonge. Merci pour toutes ces pistes à explorer.
Par exemple, il y a un résumé ici (format *.pdf, on en trouve d'autres sur le net) pour ce qui est itératif + pluri-seconds membres, mais ça risque d'être de la prise de tête pour une performance pas forcément meilleure (si tes matrices à inverser sont très grandes, N ~ plusieurs dizaines de milliers, ok mais tu risques alors d'avoir aussi un pb de stockage en RAM). Les H-matrices ne sont pas applicables partout et demandent un gros investissement en temps donc privilégie une solution sur l'étagère si tu en trouves (je n'en ai pas à te proposer là). Keep it simple, SCALAPACK je pense correctement le boulot.
Avez-vous entendu parler de Julia ? Laissez-vous tenter ...
Ce n'est plus d'actualité pour moi, mais ça m'intéresse tout de même et j'imagine que ça peut servir à d'autres aussi, donc merci beaucoup !
Inversion numérique de matrice mal conditionnée
× Après avoir cliqué sur "Répondre" vous serez invité à vous connecter pour que votre message soit publié.
× Attention, ce sujet est très ancien. Le déterrer n'est pas forcément approprié. Nous te conseillons de créer un nouveau sujet pour poser ta question.
Avez-vous entendu parler de Julia ? Laissez-vous tenter ...
Avez-vous entendu parler de Julia ? Laissez-vous tenter ...