Comment réorganiser les variables pour produire une matrice à bandes de bande passante minimale?

15

J'essaie de résoudre une équation de Poisson 2D par des différences finies. Dans le processus, j'obtiens une matrice clairsemée avec seulement variables dans chaque équation. Par exemple, si les variables étaient , alors la discrétisation donnerait: $5$ $U$

U_{i - 1, j} + U_{je + 1, j} - 4 U_{je, j} + U_{je, j - 1} + U_{je, j + 1} = F_{je, j}

$U_{i-1,j} + U_{i+1,j} -4U_{i,j} + U_{i,j-1} + U_{i,j+1} = f_{i,j}$

Je sais que je peux résoudre ce système par une méthode itérative, mais l'idée m'est venue à l'esprit que si je commandais les variables de manière appropriée, je pourrais être en mesure d'obtenir une matrice en bandes qui pourrait être résolue par une méthode directe (c.-à-d. Élimination gaussienne w / o pivotant). Est-ce possible? Existe-t-il des stratégies pour le faire pour d'autres systèmes clairsemés peut-être moins structurés?

linear-algebra pde finite-difference sparse-matrix banded-matrix Paul
la source

2

Quelque chose comme Cuthill-McKee, alors?

JM

Intéressant ... je n'ai jamais entendu parler de l'algorithme de Cuthill-McKee auparavant! :)

Paul

1

Il y a aussi un Reverse Cuthill-McKee.

Geoff Oxberry

1

J'espère qu'il ressort clairement des réponses, mais vous ne pas envie d'utiliser un solveur pour bagués ce problème, ni choisir un ordre qui réduit la bande passante. Peut-être que la question ou la réponse choisie peut être modifiée pour que cela soit clair, sinon je crains que ce mythe ne se perpétue. J'ai donné une comparaison visuelle et comparé le remplissage de scicomp.stackexchange.com/a/880/119 .

Jed Brown

@JedBrown: En fait, je ne travaille pas vraiment avec un problème de poisson, en soi ... Mon problème a une structure similaire au problème de poisson ... Les indices des variables (i et j) sont exactement les mêmes, et la matrice est dominante en diagonale avec les entrées hors diagonale (dans la même ligne) s'ajoutant exactement à la somme de l'entrée en diagonale.

Paul

13

Il s'agit d'un problème bien étudié dans le domaine des solveurs directs clairsemés. Je recommande fortement de lire l' aperçu de Joseph Liu de la méthode multifrontale afin d'avoir une meilleure idée de la façon dont les réorganisations et les supernodes affectent le temps de remplissage et de solution.

La dissection imbriquée est un moyen extrêmement courant de générer la réorganisation, et consiste essentiellement en un partitionnement récursif de graphes. Metis est la norme de facto pour le partitionnement de graphe, et vous pouvez lire quelques - unes des idées derrière elle ici . Un autre package couramment utilisé est SCOTCH , et Chaco est également important, car ses auteurs ont introduit le partitionnement graphique à plusieurs niveaux , qui est également l'idée fondamentale derrière MeTiS .

George et Liu ont montré dans leur livre classique que les solutions 2D directe rares exigent que travail et la mémoire, alors que nécessite 3d clairsemés directe travaux et mémoire. $O(n^{3/2})$ $O(n \log n)$ $O(n^2)$ $O(n^{4/3})$

Jack Poulson
la source

Avez-vous une citation pour la référence George et Liu?

Paul

Ajoutée; J'étais sur le point de sortir de la voiture lorsque je l'ai soumise pour la première fois. Je sais qu'il existe quelque part une version gratuite du livre en ligne (Jed sait où il se trouve), mais je ne l'ai pas trouvée.

Jack Poulson

J'ai mis à jour le lien pour pointer vers le PDF du livre au lieu de la critique du livre.

Jed Brown

@JedBrown C'était une excellente référence! Merci beaucoup! :)

Paul

1

@Alexander Tout le monde attribue la 3D liée à George et Liu, bien que je ne sais pas s'ils le signalent explicitement dans le livre. Cela ressort cependant de la théorie. Le séparateur de sommet minimal pour une

grille est

. La matrice dense associé à ce supernode a

entrées et nécessite

n = m \times m \times m

$n = m\times m\times m$

n^{2 / 3} = m \times m

$n^{2/3} = m\times m$

(n^{2 / 3})^{2} = n^{4 / 3}

$(n^{2/3})^2 = n^{4/3}$

(n^{2 / 3})^{3} = n^{2}

$(n^{2/3})^3 = n^2$ opérations à factoriser. Le terme logarithmique dans le cas 2D est plus subtil et est traité dans le chapitre 8 sur la dissection imbriquée, qui atteint la limite inférieure.

Jed Brown

5

Cuthill-McKee est la norme de facto pour ce que vous voulez faire. Si vous vouliez jouer avec cette méthode, il existe une implémentation facile à utiliser de l'algorithme (et son inverse) dans la Boost Graph Library (BGL), et la documentation contient des exemples d'utilisation.

Wolfgang Bangerth
la source

en fait inverser Cuhill-McKee; il donne généralement moins de remplissage. Mais un ordre de dissection imbriqué est de loin supérieur à un ordre de faible bande passante.

Arnold Neumaier

4

En parlant de méthodes multifrontales, Tim Davis , qui travaille sur les méthodes multifrontales pour la factorisation LU ( UMFPACK ), a un certain nombre de routines qui réorganiseront les matrices pour minimiser le remplissage. Vous pouvez les trouver comme ici dans le cadre de SuiteSparse . SuiteSparse utilise MeTiS.

Une autre chose à noter: dans certains problèmes, vous pouvez être intelligent pour classer les variables afin d'obtenir des modèles regroupés ou proches des modèles regroupés, ce qui peut vous éviter la peine (et le temps CPU) d'appeler ces algorithmes. Cependant, cette réorganisation intelligente nécessite un aperçu de votre part et est loin d'être aussi générale que les algorithmes de réorganisation basés sur la théorie des graphes que les gens ont mentionnés dans leurs réponses ici.

Geoff Oxberry
la source

Je vous en prie, Paul. Si vous l'aimez, votez-le.

Geoff Oxberry

3

Il y a un algorithme appelé ADI (Alternating Direction Implicit) dans les cercles de mathématiques appliquées et Split-operator dans les cercles de physique qui fait essentiellement ce que vous décrivez. C'est une méthode itérative, et elle suit cette procédure de base:

$y$ $x$
$x$ $y$
Répétez 1 et 2 jusqu'à ce que l'erreur soit aussi petite que vous le souhaitez.

Je ne connais pas la complexité formelle de cet algorithme, mais je l'ai trouvé converger en moins d'itérations que des choses comme Jacobi et Gauss-Seidel chaque fois que je l'ai utilisé.

Dan
la source

2

Si vous décidez de suivre la voie de fractionnement de l'opérateur, vous devrez faire attention à ce que le fractionnement de l'opérateur soit connu pour entraîner des erreurs dans les solutions en régime permanent dans certains cas. (Un de mes camarades de laboratoire a développé un moyen de surmonter cette difficulté, mais je ne pense pas qu'il l'ait encore publié.) De plus, le fractionnement d'opérateur est connu pour provoquer des erreurs numériques. Il existe des moyens bien établis d'estimer ces erreurs a posteriori ; Don Estep a fait un excellent travail dans ce domaine.

Geoff Oxberry

@GeoffOxberry On dirait que vous faites référence à un fractionnement différent. Vous pouvez utiliser ADI dans un schéma entièrement implicite qui n'a pas d'erreur de fractionnement car il résout réellement le système. Il existe également des méthodes IMEX qui contrôlent rigoureusement les erreurs de fractionnement.

Jed Brown

x

$x$

y

$y$

Je n'ai jamais entendu parler de la séparation de Godunov et Strang. J'ai tendance à diviser mon opérateur avec la formule Baker-Campbell-Hausdorf. Est-ce la même chose?

Dan

Comment réorganiser les variables pour produire une matrice à bandes de bande passante minimale?

Réponses: