Meilleur choix de solveur pour un grand système symétrique clairsemé (mais pas défini positif)

Je travaille actuellement sur la résolution de très grands systèmes symétriques (mais pas définis positifs), générés par certains algorithmes. Ces matrices ont une belle densité de blocs qui peut être utilisée pour la résolution parallèle. Mais je ne peux pas décider si je dois utiliser une approche directe (comme Multi-frontale) ou itérative (GMRES ou MINRES préconditionnée). Toutes mes études montrent que les solveurs itératifs (même avec une convergence assez rapide de 7 itérations internes) ne parviennent pas à battre l'opérateur '\' direct dans MATLAB. Mais en théorie, les méthodes directes sont censées être plus coûteuses. Comment cela se passe-t-il? Existe-t-il un document ou un papier à jour pour un tel cas? Puis-je utiliser la parcimonie par blocs dans des systèmes parallèles en utilisant des méthodes directes aussi efficacement que des solveurs itératifs flexibles comme GMRES?

Le solveur direct clairsemé MUMPS peut gérer des systèmes indéfinis symétriques et est disponible gratuitement ( http://graal.ens-lyon.fr/MUMPS/ ). Ian Duff était l'un des auteurs de MUMPS et de MA57, les algorithmes présentent donc de nombreuses similitudes.

MUMPS a été conçu pour les ordinateurs parallèles à mémoire distribuée, mais il fonctionne également bien sur les machines à processeur unique. Si vous le liez à une bibliothèque BLAS multithread, vous pouvez obtenir des accélérations raisonnables sur la mémoire partagée, les processeurs multicœurs.

Vous n'avez pas dit quelle taille est «très grande», mais MUMPS a également une capacité hors cœur pour gérer les problèmes où la matrice factorisée ne tient pas dans la mémoire disponible.

Le problème général dont souffrent les solveurs directs est le phénomène de remplissage, ce qui signifie que l'inverse d'une matrice clairsemée peut être dense. Cela entraîne d'énormes besoins en mémoire si la structure de la matrice n'est pas "adaptée".

Il existe des tentatives pour contourner ces problèmes, et la fonction par défaut de luMATLAB en emploie quelques-uns. Voir http://www.mathworks.de/de/help/matlab/ref/lu.html pour un aperçu des permutations, de la mise à l'échelle diagonale et autres. Il en va de même bien sûr pour les packages plus avancés comme MUMPS ( http://graal.ens-lyon.fr/MUMPS/ ).

En général, cependant, si votre problème est suffisamment important, vous atteindrez cette limite de mémoire et vous devrez utiliser des méthodes itératives (Krylov).

Si votre problème est symétrique et indéfini, MINRES peut être le choix évident. Notez, cependant, que GMRES et MINRES font la même chose en arithmétique exacte si le problème est symétrique, mais GMRES a tendance à souffrir moins de la perte d'orthogonalité. Par conséquent, certains considèrent GMRES comme "la meilleure mise en œuvre de MINRES".

Dans tous les cas, vous profiterez probablement du préconditionnement de votre système. Si vous ne voulez pas consacrer de temps à la personnalisation d'un préconditionneur, un préconditionneur de factorisation LU (ILU) incomplet peut déjà vous mener quelque part.

Tout solveur itératif ne peut battre les méthodes directes que si le problème est suffisamment important (grand, dépend de plusieurs facteurs tels que le stockage requis, l'efficacité de la mise en œuvre). Et toutes les méthodes krylov (par exemple GMRES) ne sont bonnes que si vous utilisez un préconditionnement approprié (en pratique). Si vous n'utilisez aucun préconditionnement, les méthodes krylov ne sont généralement pas utiles. Moi aussi, je travaille avec des matrices clairsemées de blocs (celles-ci proviennent d'applications PDE) et j'ai observé que le solveur de krylov préconditionné par bloc (chevauchement d'additifs schwarz) (soit redémarré GMRES ou BiCG-Stab) couplé avec une correction de grille grossière (ou une multigrille) est évolutif pour ces type de matrices.

L'opérateur Matlab '\' est très efficace grâce à une programmation de premier ordre. Vous pouvez voir une partie de l'idée et comment ils ont utilisé tous les raccourcis possibles dans le livre de Timothy Davis.

Dans matlab, vous pouvez utiliser gmres, qui est bien développé et stable. Les mineurs, qui devraient théoriquement être idéaux pour votre cas, peuvent ne pas être aussi fiables (du moins mon expérience le dit). Vous devriez obtenir une efficacité égale ou supérieure de matlab gmres si

1. Votre système est suffisamment grand (au moins quelques milliers par quelques milliers).
2. Si vous utilisez le bon type de paramètres (RESTART, MAXIT, X0). Aucune directive claire n'est disponible pour cela. Utilisez votre expérience.
3. Utilisez un bon préconditionneur. Vous pouvez utiliser ILU ou encore moins cher, un bloc Jacobi. Cela réduira considérablement l'effort.
4. LE PLUS IMPORTANT: Si votre matrice est clairsemée, utilisez le format matlab clairsemé. Matlab gmres est idéalement construit pour cela. Cela réduira les coûts dans une large mesure.

Pour des systèmes encore plus grands, utilisez un outil comme PETSc.

Si votre matrice a une dimension au milieu des dizaines de milliers ou moins, utilisez une méthode directe, bien qu'il n'y ait pas beaucoup de méthodes directes disponibles gratuitement pour les systèmes symétriques indéfinis (en fait, aucune que je sache n'est open source). Il y a MA57 de la HSL, mais ce n'est gratuit que pour une utilisation académique. Vous pouvez certainement ignorer la symétrie et utiliser UMFPACK .

Aux alentours de quelques dizaines de centaines de dimensions, l'utilisation de la mémoire d'une méthode directe commence à dépasser ce qu'un ordinateur de bureau typique peut raisonnablement gérer, donc à moins d'avoir une machine à mémoire partagée robuste, vous devrez passer à des méthodes itératives. Pour les problèmes indéfinis, vous pouvez spécialiser le BiCG (gradient biconjugué) pour les systèmes symétriques, bien qu'une panne soit possible. Il existe un MINRES-QLP récemment publié pour les systèmes symétriques qui offre une plus grande stabilité numérique.

Les deux méthodes nécessitent des implémentations plutôt différentes, car pour les méthodes directes, vous devez réellement former la matrice, tandis que dans l'approche itérative, vous ne formerez généralement pas la matrice explicitement.

Il existe un certain nombre de raisons pour lesquelles une approche peut être plus rapide que l'autre, en particulier en fonction de la dimension de la matrice. Pour les systèmes très mal conditionnés, les méthodes itératives peuvent stagner assez mal. Pour les matrices moins clairsemées, les méthodes directes finissent par créer beaucoup de remplissage, ce qui ralentit beaucoup les choses. De plus, les méthodes directes dans Matlab sont hautement optimisées (il utilise UMFPACK ou MA57 en interne), tandis que les méthodes itératives sont généralement codées directement dans Matlab, et il y a moins de possibilités d'exploiter BLAS de niveau 3 car les goulots d'étranglement sont l'application de matvecs et de produits scalaires.