Inversion de matrice en grande dimension

Bonsoir,

Je recherche un algorithme déjà existant ou une méthode à coder pour inverser une matrice de grande dimension (8000 par 8000). J'ai déjà des solutions en $O(n^3)$, mais cela fait dans les 3 jours de temps de calcul (avec ce que j'ai déjà vu), et je ne suis pas sûr d'avoir une garantie concernant le résultat. Dans l'absolu le fait de consacrer 3 jours de temps CPU à cette opération ne me dérange pas (les week-end sont faits pour ça) , mais j'aimerais évidemment être sûr du résultat final.

Quelqu'un aurait-il déjà eu ce problème et rencontré des méthodes efficaces pour réaliser cette inversion ?
Je précise, à tout hasard pour mes amis numériciens, que la matrice que je souhaite inverser est symétrique, au cas où cette information pourrait améliorer les choses.

Amicalement,

Réponses

  • Serait-elle creuse par miracle ? (on ne sait jamais)
  • Malheureusement non, faut pas rêver ... ;)
  • Tu ne fera jamais mieux que $O(n^3)$ en utilisant l' algorithme de Gauss par exemple
  • tu peux la preconditionner par une matrice creuse?
  • Bonjour,

    Quelques questions à se poser (en plus de celle du barbu rasé):
    - s'agit-il vraiment d'inverser la matrice A ou bien de résoudre un système du type AX = X0 ?
    - y a t-il, a priori, des problèmes de stabilité ?
    - si, en plus d'être symétrique, la matrice est définie positive, alors la méthode de Choleski est particulièrement indiquée.

    L'origine physique du problème peut suggérer des éléments de réponses.

    Cordialement
  • Pilz : je sais que Gauss est en $O(n!)$ et que donc ce n'est absolument pas ce que je dois utiliser (je n'ose même pas utiliser Stirling pour voir à peu près combien vaut $8000!$ ...). Je cherche des méthodes "nouvelles et efficaces" si je puis dire, pour peu qu'elles existent.

    Gecko : préconditionner je veux bien, mais qu'est-ce que ça pourrait apporter du point de vue temps de calcul de l'inverse ? Cette question est naïve et non pas agressive.
  • GPP29 :
    - Il s'agit vraiment de l'inverser
    - Qu'entends-tu par problème de stabilité ?
    - La matrice est bien définie positive. Mais Choleski en grande dimension ?

    Enfin, l'origine mathématique plutôt et non pas physique de ce problème est une régression linéaire à 8000 régresseurs.
  • Vous programmez en quoi ? La routine Lapack est ce qui se fait de mieux en ce moment pour tout ce qui est calcul matriciel, elle peut être utilisée sous Fortran, C et elle est intégrée à Matlab.

    La méthode de Gauss est bien en O(n^3) et non en O(n!) comme tu le dis. Cholesky est plus stable que Gauss mais c'est aussi du O(n^3).

    Pour une régression linéaire je crois qu'on se ramène à la résolution d'un système et non pas l'inversion d'une matrice, ce qui permet de gagner pas mal de temps. Une factorisation QR (toujours en O(n^3) ) peut aussi permetre d'augmenter la stabilité.
  • Pour Kuja,

    Problème de stabilité : pour faire simple et intuitif, on parle de problème de stabilité lorsque de petites variations des données initiales produisent de grandes variations des résultats. Dans un tel cas, l'accumulation des erreurs d'arrondis au cours des calculs peut produire des résultats dénués de sens. C'est d'autant plus à craindre que les matrices sont grandes.

    Concernant Cholesky en grande dimension : des études menées à l'Université de Stuttgart (Intitut für Statik und Dynamic) ont comparé l'efficacité de divers algorithmes, sur des problèmes issus de la mécanique. Les méthodes de type Cholesky, avec différentes stratégies de stockage des données, se sont avérées de très loin les meilleures. Elles conviennent particulièrement bien aux grandes dimensions et sont remarquables de stabilité. Avec les éléments finis, on atteint couramment des centaines de milliers de variables.
  • Par le préconditionnement on est en O(n^3) et l'inverse est approché (de toute façon c'est toujours le cas) mais d'après ce que je viens de lire...
    Ce n'est pas le bon choix
    Bon courage

    Que le CPU chauffe bien!
  • O(n^4) pour une inversion...
    O(n^3) pour une solution du système
  • Bonjour,

    Merci à tous pour vos réponses.

    En ce qui concerne Gauss, j'étais persuadé que c'était en $O(n!)$, autant pour moi.
    D'autre part, autre énorme bourde de ma part : comme le faisait remarquer seb, évidemment en régression linéaire on ne demande que de réussir à résoudre un système linéaire.

    Donc je résume :

    - J'ai en fait un système linéaire de grande taille (8000 par 8000) à résoudre
    - J'aimerais avoir la méthode la plus stable
    - A priori, il faut que j'aille jeter un oeil à Choleski

    Quelqu'un aurait-il une référence pour Cholesky en grande dimension ?
    Le fait qu'il s'agisse d'un système linéaire et non pas d'inversion pure permet il d'envisager une stratégie plus efficace ?

    Merci encore et désolé pour mes bourdes.
  • Sinon, avec l'algorithme de Strassen, on peut inverser de manière exacte une matrice en $\mathcal{O}(n^{\alpha})$ opérations, avec $\alpha = \frac{ln(7)}{ln(2)} \approx 2.81...$

    Ou alors, on a aussi des méthodes itératives (Jacobi, Gauss-Seidel,...) qui vont plus vite. Le problème, c'est qu'elles ne donnent pas l'inverse de la matrice de manière exacte, mais fournissent juste des suites qui convergent (vite) vers l'inverse de la matrice. Alors si tu n'as besoin que d'un résultat numérique à $10^{-3}$ près, par exemple, c'est certainement le meilleur moyen.

    Pour plus de détails, tu peux aller voir dans "Matrices", de Denis Serre.
  • > kuja, je sais que ca n'a rien a voir, mais pas besoin de stirling pour vir combien font 8000! : c'est une bonne adresse a connaitre quand on a oublie sa calculette : <a href=" http://www.swox.com/gmp/#TRY"&gt; http://www.swox.com/gmp/#TRY</a&gt; permet d'essayer en ligne la bibliotheque C GMP, et donc de faire de "petits" calcul, ce qui s'aver parfois bien pratique !! en l'occurence ca me sort 8000! en 8 ms, et ca fait plutot gros...<BR><BR><BR>
  • Guego : Pour Strassen, tu l'as déjà programmé et fait marcher pour de vrai (i.e pas sur une matrice 4*4) ? Parce que c'est typiquement le genre de résutat théorique agréable qui est abominable à implémenter, du moins pour ce que j'en ai vu.
  • Oui, je l'ai déjà fait, et c'est pas si immonde que ça. En tout cas, beaucoup moins que certains autres algos d'analyse numérique.
  • Merci Guego, je vais y jeter un oeil et reviendrai poser des questions si besoin est.
  • Juste par curiosité, dans quel langage l'avais tu codé ? Et sais tu si ce type d'algo est implémenté dans les standards d'analyse numérique (Matlab, Scilab ...) ?
  • Je vois parler de Gauss. Règle de base de l'analyse numérique matricielle : on oublie Gauss, méthode instable. Essayez donc de faire du Gauss avec une matrice de Hilbert ...
  • Eric : C'est bien pour ça que je parlais de la méthode QR

    Je reste quand même pour l'utilisation de la librairie Lapack qui a été largement optimisée.
  • Strassen c'est pas trop dur à programmer je l' ai fait en Maple. Cependant mais si c'est théoriquement en $O(n^2.81)$ ce n'est absolument pas plus rapide que les méthodes d' inversion usuelles tellement la constante cachée dans le $O$ est importante.

    A noter que l' on ne connais pas la complexité exacte de l' inversion d 'une matrice , on sait montrer que ce cout est équivalent au cout de multiplication de 2 matrices...
    Le meilleur algorithme trouvé à ce jour est en $O(n^2.30)$
  • Pour mes complexité il faut lire $O(n^{2.81})$ et $O(n^{2.30})$
  • Bonsoir,

    juste un petit mot pour remercier tous les intervenants du topic, car j'ai réussi à résoudre mon problème.
    Un coup de Choleski au départ (une vingtaine de minute sur ma machine pour une matrice 8000 par 8000) avec les routines Lapack, puis la résolution de deux systèmes triangulaires (une dizaine de minutes chacun) et le tour est joué, avec une erreur en $10^{-13}$ sur les résultats. Bref, le bonheur pour moi (enfin plutôt pour ma thèse).

    Un grand merci à tous, et je recommande vivement Cholesky à tout utilisateur qui rencontrerait le même problème que moi, ça marche du tonnerre.

    Amicalement,
Connectez-vous ou Inscrivez-vous pour répondre.