Comment l'opérateur de barre oblique inversée MATLAB résout-il

36

Je comparais quelques-uns de mes codes avec les codes MATLAB "en stock". Je suis surpris des résultats.

J'ai couru un exemple de code (Sparse Matrix)

n = 5000;
a = diag(rand(n,1));
b = rand(n,1);
disp('For a\b');
tic;a\b;toc;
disp('For LU');
tic;LULU;toc;
disp('For Conj Grad');
tic;conjgrad(a,b,1e-8);toc;
disp('Inv(A)*B');
tic;inv(a)*b;toc;

Résultats :

    For a\b
    Elapsed time is 0.052838 seconds.

    For LU
    Elapsed time is 7.441331 seconds.

    For Conj Grad
    Elapsed time is 3.819182 seconds.

    Inv(A)*B
    Elapsed time is 38.511110 seconds.

Pour la matrice dense:

n = 2000;
a = rand(n,n);
b = rand(n,1);
disp('For a\b');
tic;a\b;toc;
disp('For LU');
tic;LULU;toc;
disp('For Conj Grad');
tic;conjgrad(a,b,1e-8);toc;
disp('For INV(A)*B');
tic;inv(a)*b;toc;

Résultats:

For a\b
Elapsed time is 0.575926 seconds.

For LU
Elapsed time is 0.654287 seconds.

For Conj Grad
Elapsed time is 9.875896 seconds.

Inv(A)*B
Elapsed time is 1.648074 seconds.

Comment diable est un \ b si génial?

Enquête
la source
1
La barre oblique inverse intégrée de MATLAB, autrement dit le résolveur direct pour un système d'équations linéaires, utilise la méthode Multifrontal pour les matrices creuses, c'est pourquoi A \ B est si génial.
Shuhao Cao
1
Il utilise le code de Tim Davis disponible à l' adresse cise.ufl.edu/research/sparse . Aussi la génialité disparaît quand vous avez un problème non trivial
Stali
1
Qu'est ce que "LULU"? Pourquoi pensez-vous que c'est une bonne implémentation d'une factorisation LU et une résolution directe ultérieure?
Jed Brown
3
@Nunoxic Quelle implémentation? L'avez-vous écrit vous-même? L'algèbre linéaire dense à haute performance, bien que généralement bien comprise sur le plan algorithmique, n'est pas facile à mettre en œuvre efficacement sur du matériel moderne. Les meilleures implémentations de BLAS / Lapack devraient être proches du maximum pour une matrice de cette taille. De plus, à partir de vos commentaires, j'ai l’impression que vous pensez que LU et Gaussian Elimination sont des algorithmes différents.
Jed Brown
1
Il appelle un code Fortran écrit avec Intel MKL.
Enquête le

Réponses:

37

Dans Matlab, la commande '\' appelle un algorithme qui dépend de la structure de la matrice A et comprend des vérifications (faible temps système) sur les propriétés de A.

  1. Si A est rare et en bandes, utilisez un solveur en bandes.
  2. Si A est une matrice triangulaire supérieure ou inférieure, utilisez un algorithme de substitution en arrière.
  3. Si A est symétrique et que les éléments diagonaux positifs sont réels, essayez une factorisation de Cholesky. Si A est rare, utilisez d'abord le réordonnancement pour minimiser le remplissage.
  4. Si aucun des critères ci-dessus n’est rempli, effectuez une factorisation triangulaire générale en utilisant l’élimination gaussienne avec pivot partiel.
  5. Si A est rare, utilisez la bibliothèque UMFPACK.
  6. Si A n'est pas carré, utilisez des algorithmes basés sur la factorisation QR pour des systèmes indéterminés.

Pour réduire les frais généraux, il est possible d'utiliser la commande linsolve dans Matlab et de sélectionner vous-même un solveur approprié parmi ces options.

Allan P. Engsig-Karup
la source
En supposant que je traite avec une matrice dense non structurée 10000x10000 avec tous les éléments non nuls (niveau de densité élevé), quel serait mon meilleur choix? Je veux isoler cet algorithme 1 qui fonctionne pour les matrices denses. Est-ce LU, QR ou Gaussian Elimination?
Enquête le
1
Cela ressemble à une étape 4 où l’élimination gaussienne est invoquée, ce qui correspond au cas le plus général dans lequel aucune structure de A ne peut être exploitée pour améliorer les performances. Il s’agit donc d’une factorisation de la LU et de la suivante, suivie d’une étape de substitution en arrière.
Allan P. Engsig-Karup
Merci! Je pense que cela me donne une direction à penser. À l’heure actuelle, l’élimination gaussienne est la meilleure solution pour résoudre de tels problèmes non structurés, est-ce exact?
Enquête le
37

Si vous voulez voir ce qui a\bfait pour votre matrice particulière, vous pouvez définir spparms('spumoni',1)et déterminer exactement quel algorithme vous a impressionné. Par exemple:

spparms('spumoni',1);
A = delsq(numgrid('B',256));
b = rand(size(A,2),1);
mldivide(A,b);  % another way to write A\b

va sortir

sp\: bandwidth = 254+1+254.
sp\: is A diagonal? no.
sp\: is band density (0.01) > bandden (0.50) to try banded solver? no.
sp\: is A triangular? no.
sp\: is A morally triangular? no.
sp\: is A a candidate for Cholesky (symmetric, real positive diagonal)? yes.
sp\: is CHOLMOD's symbolic Cholesky factorization (with automatic reordering) successful? yes.
sp\: is CHOLMOD's numeric Cholesky factorization successful? yes.
sp\: is CHOLMOD's triangular solve successful? yes.

je peux donc voir que "\" a fini par utiliser "CHOLMOD" dans ce cas.

dranxo
la source
3
+1 pour les nouveaux paramètres MATLAB dont je n'avais jamais entendu parler. Bien joué, monsieur.
Geoff Oxberry
2
Hey, merci! C'est dans help mldivide.
dranxo
16

Pour les matrices creuses, Matlab utilise UMFPACK pour l' \opération " ", qui, dans votre exemple, parcourt les valeurs de a, les inverse et les multiplie par les valeurs de b. Pour cet exemple, vous devriez utiliser b./diag(a), ce qui est beaucoup plus rapide.

Pour les systèmes denses, l'opérateur de barre oblique inverse est un peu plus compliqué. Une brève description de ce qui est fait quand est donnée ici . Selon cette description, dans votre exemple, Matlab résoudrait l’ a\butilisation de la substitution arrière. Pour les matrices carrées générales, une décomposition de LU est utilisée.

Pedro
la source
Regd. Sparsity, l’inv d’une matrice de diagnostic serait simplement l’inverse des éléments diagonaux pour que b./diag(a) fonctionne, mais a \ b fonctionne également pour les matrices rares. Pourquoi linsolve ou LULU (ma version optimisée de LU) n’est-il pas plus rapide que a \ b pour les matrices denses dans ce cas.
Enquête le
@Nunoxic Votre LULU fait-il quelque chose pour détecter la diagonalité ou la triangularité de la matrice dense? Sinon, cela prendra aussi longtemps avec chaque matrice, quel que soit son contenu ou sa structure.
Pedro
Quelque peu. Mais, avec les drapeaux OPT de linsolve, j’ai tout défini pour définir la structure. Pourtant, a \ b est plus rapide.
Enquête le
@Nunoxic, chaque fois que vous appelez une fonction utilisateur, vous créez des frais généraux. Matlab fait tout pour la barre oblique inverse en interne, par exemple la décomposition puis l’application ultérieure du côté droit, avec très peu de frais généraux, et sera donc plus rapide. De plus, dans vos tests, vous devez utiliser plus d’un appel pour obtenir des horaires fiables, par exemple tic; for k=1:100, a\b; end; toc.
Pedro
5

En règle générale, si vous avez une matrice creuse de complexité raisonnable (c’est-à-dire que ce n’est pas nécessairement le gabarit à 5 points, mais peut en fait être une discrétisation des équations de Stokes pour laquelle le nombre de nonzeros par ligne est beaucoup plus grand que 5), un solveur direct maigre, tel que UMFPACK, bat généralement un solveur de Krylov itératif si le problème ne dépasse pas environ 100 000 inconnus.

En d'autres termes, pour la plupart des matrices creuses résultant de la discrétisation 2D, un résolveur direct est la solution la plus rapide. Il est impératif d'utiliser un résolveur itératif uniquement pour les problèmes 3D où vous obtenez rapidement plus de 100 000 inconnus.

Wolfgang Bangerth
la source
3
Je ne vois pas comment cela répond à la question, mais je conteste également le principe. Il est vrai que les solveurs directs ont tendance à bien fonctionner pour les problèmes 2D de taille modeste, mais les solveurs directs ont tendance à souffrir en 3D bien avant 100 000 inconnues (les séparateurs de vertex sont beaucoup plus grands qu'en 2D). De plus, j'affirme que dans la plupart des cas (opérateurs elliptiques, par exemple), les solveurs itératifs peuvent être utilisés plus rapidement que les solveurs directs, même pour des problèmes 2D de taille moyenne, mais cela peut demander un effort considérable (par exemple, utiliser les bons ingrédients pour que multigrid fonctionne) .
Jed Brown
1
Si vos problèmes sont raisonnablement modestes et que vous ne souhaitez pas concevoir de solutionneurs implicites, ou si vos problèmes sont très difficiles (comme Maxwell à haute fréquence) et si vous ne souhaitez pas consacrer votre carrière à la conception de solutions correctives, convenez que les solveurs directs épars sont un excellent choix.
Jed Brown
En fait, mon problème est de concevoir un bon résolveur dense. Je n'ai pas d'application particulière en tant que telle (donc, pas de structure). Je voulais modifier quelques-uns de mes codes et les rendre concurrentiels par rapport à ce qui est actuellement utilisé. J'étais simplement curieux de savoir sur a \ b et sur le fait que cela compense toujours la merde de mon code.
Enquête le
@JedBrown: Oui, peut-être qu'avec beaucoup d'efforts, on peut obtenir des solveurs itératifs pour battre un solveur direct, même pour des problèmes mineurs en 2D. Mais pourquoi le faire? Pour les problèmes avec moins de 100 000 inconnus, les solveurs directs maigres en 2D sont assez rapides. Ce que je voulais dire, c'est: pour de si petits problèmes, ne passez pas votre temps à bricoler et à trouver la meilleure combinaison de paramètres - prenez simplement la boîte noire.
Wolfgang Bangerth
Il existe déjà une différence d’ordre de grandeur entre un multigrille géométrique de Cholesky clairsemé et un multigrille géométrique (à base de matrice) pour des problèmes 2D "faciles" avec 100 000 dof utilisant un gabarit à 5 points (environ 0,5 seconde contre environ 0,05 seconde). Si votre gabarit utilise des voisins proches (par exemple, discrétisation du quatrième ordre, Newton pour certains choix de rhéologie non linéaire, de stabilisation, etc.), la taille des séparateurs de vertex double environ, de sorte que le coût de la résolution directe est multiplié par 8 (le coût est plus élevé). dépendant du problème pour MG). Avec beaucoup de pas de temps ou de boucles d'optimisation / UQ, ces différences sont significatives.
Jed Brown