Foreplay

La solution numérique des systèmes d'équations linéaires est une étape indispensable dans de nombreux domaines des mathématiques appliquées, de l'ingénierie et de l'industrie informatique, qu'il s'agisse de travailler avec des graphiques, de calculer l'aérodynamique d'un avion ou d'optimiser la logistique. La «machine» désormais à la mode n'aurait pas non plus beaucoup progressé sans elle. De plus, la solution des systèmes linéaires, en règle générale, consomme le plus grand pourcentage de tous les coûts informatiques. Il est clair que ce composant critique nécessite une optimisation de la vitesse maximale.

Travaille souvent avec les soi-disant matrices clairsemées - celles dont les zéros sont des ordres de grandeur supérieurs aux autres valeurs. Ceci, par exemple, est inévitable si vous avez affaire à des équations aux dérivées partielles ou à d'autres processus dans lesquels les éléments émergents dans les relations linéaires de définition ne sont connectés qu'avec des «voisins». Voici un exemple possible d'une matrice clairsemée pour l'équation de Poisson unidimensionnelle connue en physique classique $$$- \ phi ^ {''} = f$ sur une grille uniforme (oui, bien qu'il n'y ait pas tant de zéros dedans, mais lors du broyage de la grille, ils seront sains!):

$$

$$\ begin {pmatrix} 2 & -1 & 0 & 0 & 0 \\ -1 & 2 & -1 & 0 & 0 \\ 0 & -1 & 2 & -1 & 0 \\ 0 & 0 & -1 & 2 & -1 \\ 0 & 0 & 0 & -1 & 2 \ end {pmatrix}$$

Les matrices qui leur sont opposées - celles dans lesquelles elles ne font pas attention au nombre de zéros et prennent en compte tous les composants sans exception - sont appelées denses.

Les matrices éparses sont souvent, pour diverses raisons, présentées dans un format de colonne compressée - Colonne creuse compressée (CSC). Cette description utilise deux tableaux entiers et un point flottant. Laissez la matrice tout avoir $$$nnz_A$ non nul et $$$N$ colonnes. Les éléments de la matrice sont répertoriés dans des colonnes de gauche à droite, sans exception. Premier tableau $$$i_A$ longueurs $$$nnz_A$ contient des numéros de ligne de composants de matrice non nuls. Deuxième tableau $$$j_A$ longueurs $$$N + 1$ contient ... non, pas des numéros de colonnes, car alors la matrice serait écrite au format de coordonnées (coordonnées) ou ternaire (triplet). Et le deuxième tableau contient les numéros de série des composants de matrice avec lesquels les colonnes commencent, y compris une colonne factice supplémentaire à la fin. Enfin, le troisième tableau $$$v_A$ longueurs $$$nnz_A$ contient déjà les composants eux-mêmes, dans l'ordre correspondant au premier tableau. Ici, par exemple, sous l'hypothèse que la numérotation des lignes et des colonnes de matrices commence à zéro, pour une matrice particulière

$$

$$A = \ begin {pmatrix} 0 & 1 & 0 & 4 & -1 \\ 7 & 0 & 2.1 & 0 & 3 \\ 0 & 0 & 0 & 10 & 0 \ end {pmatrix}$$

les tableaux seront $$$i_A = \ {1, 0, 1, 0, 2, 0, 1 \}$ , $$$j_A = \ {0, 1, 2, 3, 5, 7 \}$ , $$$v_A = \ {7, 1, 2.1, 4, 10, -1, 3 \}$ .

Les méthodes de résolution de systèmes linéaires sont divisées en deux grandes classes - directes et itératives. Les lignes droites sont basées sur la possibilité de représenter la matrice comme un produit de deux matrices plus simples, afin de diviser ensuite la solution en deux étapes simples. Les itératifs utilisent les propriétés uniques des espaces linéaires et travaillent sur l'ancien comme méthode mondiale d'approximation successive des inconnues à la solution souhaitée, et dans le processus de convergence lui-même, les matrices ne sont généralement utilisées que pour les multiplier par des vecteurs. Les méthodes itératives sont beaucoup moins chères que les méthodes directes, mais fonctionnent lentement sur des matrices mal conditionnées. Lorsque la fiabilité du béton armé est importante - utilisez des lignes droites. Me voici ici et je veux aborder un peu.

Disons que pour une matrice carrée $$$A$ la décomposition du formulaire est à notre disposition $$$A = LU$ où $$$L$ et $$$U$ - respectivement, les matrices triangulaires inférieure et triangulaire supérieure, respectivement. Le premier signifie qu'il a un zéro au-dessus de la diagonale, le second signifie qu'il est en dessous de la diagonale. Comment exactement nous avons obtenu cette décomposition - nous ne sommes pas intéressés maintenant. Voici un exemple de décomposition simple:

$$

$$\ begin {pmatrix} 1 & -1 & -1 \\ 2 & - 1 & -0.5 \\ 4 & -2 & -1.5 \ end {pmatrix} = \ begin {pmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ 4 & 2 & 1 \ end {pmatrix} \ begin {pmatrix} 1 & -1 & -1 \\ 0 & 1 & 1.5 \\ 0 & 0 & -0.5 \ end {pmatrix}$$

Comment dans ce cas résoudre le système d'équations $$$Ax = f$ par exemple, avec le côté droit $$$f = \ begin {pmatrix} 4 \\ 2 \\ 3 \ end {pmatrix}$ ? La première étape est un mouvement direct (résolution avant = substitution avant). Nous dénotons $$$y: = Ux$ et travailler avec le système $$$Ly = f$ . Parce que $$$L$ - triangulaire inférieur, successivement dans le cycle on retrouve tous les composants $$$y$ de haut en bas:

$$

$$\ begin {pmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ 4 & 2 & 1 \ end {pmatrix} \ begin {pmatrix} y_1 \\ y_2 \\ y_3 \ end {pmatrix} = \ begin {pmatrix} 4 \\ 2 \\ 3 \ end {pmatrix} \ implique y = \ begin {pmatrix} 4 \\ -6 \\ -1 \ end {pmatrix}$$

L'idée centrale est que, en trouvant $$$i$ composants vectoriels $$$y$ il est multiplié par une colonne avec le même numéro de matrice $$$L$ qui est ensuite soustrait du côté droit. La matrice elle-même semble s'effondrer de gauche à droite, diminuant de taille à mesure que l'on trouve de plus en plus de composants vectoriels $$$y$ . Ce processus est appelé élimination des colonnes.

La deuxième étape est un mouvement en arrière (résolution en arrière = substitution en arrière). Avoir trouvé un vecteur $$$y$ nous décidons $$$Ux = y$ . Ici, nous parcourons déjà les composants de bas en haut, mais l'idée reste la même: $$$i$ e colonne est multipliée par le composant qui vient d'être trouvé $$$x_i$ et se déplace vers la droite, et la matrice s'effondre de droite à gauche. L'ensemble du processus est illustré pour la matrice mentionnée dans l'exemple de l'image ci-dessous:

$$

$$\ small \ begin {pmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ 4 & 2 & 1 \ end {pmatrix} \ begin {pmatrix} y_1 \\ y_2 \\ y_3 \ end {pmatrix} = \ begin {pmatrix} 4 \\ 2 \\ 3 \ end {pmatrix} \ implique \ begin {pmatrix} 1 & 0 \\ 2 & 1 \ end {pmatrix} \ begin {pmatrix} y_2 \\ y_3 \ end {pmatrix } = \ begin {pmatrix} -6 \\ -13 \ end {pmatrix} \ implique \ begin {pmatrix} 1 \ end {pmatrix} \ begin {pmatrix} y_3 \ end {pmatrix} = \ begin {pmatrix} -1 \ end {pmatrix}$$

$$

$$\ small \ begin {pmatrix} 1 & -1 & -1 \\ 0 & 1 & 1.5 \\ 0 & 0 & -0.5 \ end {pmatrix} \ begin {pmatrix} x_1 \\ x_2 \\ x_3 \ end { pmatrix} = \ begin {pmatrix} 4 \\ -6 \\ -1 \ end {pmatrix} \ Rightarrow \ begin {pmatrix} 1 & -1 \\ 0 & 1 \ end {pmatrix} \ begin {pmatrix} x_1 \ \ x_2 \ end {pmatrix} = \ begin {pmatrix} 6 \\ -9 \ end {pmatrix} \ implique \ begin {pmatrix} 1 \ end {pmatrix} \ begin {pmatrix} x_1 \ end {pmatrix} = \ begin {pmatrix} -3 \ end {pmatrix}$$

et notre décision sera $$$x = \ begin {pmatrix} -3 \\ - 9 \\ 2 \ end {pmatrix}$ .

Si la matrice est dense, c'est-à-dire qu'elle est complètement représentée sous la forme d'une sorte de tableau, unidimensionnel ou bidimensionnel, et l'accès à un élément spécifique en elle a lieu au fil du temps $$$O (1)$ , alors une procédure de solution similaire avec la décomposition existante n'est pas difficile et facile à coder, nous n'y perdrons donc même pas de temps. Et si la matrice est clairsemée? Ici, en principe, n'est pas non plus difficile. Voici le code C ++ pour le mouvement vers l'avant dans lequel la solution $$$x$ écrit pour l'endroit sur le côté droit, sans vérifier l'exactitude des données d'entrée (les tableaux CSC correspondent à la matrice triangulaire inférieure):

Algorithme 1:

void forward_solve (const std::vector<size_t>& iA, const std::vector<size_t>& jA, const std::vector<double>& vA, std::vector<double>& x) { size_t j, p, n = x.size(); for (j = 0; j < n; j++) //    { x[j] /= vA[jA[j]]; for (p = jA[j]+1; p < jA[j+1]; p++) x[iA[p]] -= vA[p] * x[j] ; } } 

Toutes les discussions ultérieures ne concerneront que la course avant pour résoudre le système triangulaire inférieur $$$Lx = f$ .

Cravate

Et si le côté droit, c'est-à-dire vecteur à droite du signe égal $$$Lx = f$ At-il un grand nombre de zéros? Ensuite, il est logique d'ignorer les calculs associés aux positions nulles. Les modifications du code dans ce cas sont minimes:

Algorithme 2:

 void fake_sparse_forward_solve (const std::vector<size_t>& iA, const std::vector<size_t>& jA, const std::vector<double>& vA, std::vector<double>& x) { size_t j, p, n = x.size(); for (j = 0; j < n; j++) //    { if(x[j]) { x[j] /= vA[jA [j]]; for (p = jA[j]+1; p < jA[j+1]; p++) x[iA[p]] -= vA[p] * x[j]; } } } 

La seule chose que nous avons ajoutée est l' if , dont le but est de réduire le nombre d'opérations arithmétiques à leur nombre réel. Si, par exemple, tout le côté droit est constitué de zéros, alors vous n'avez pas besoin de compter quoi que ce soit: la solution sera le côté droit.

Tout a l'air génial et bien sûr, cela fonctionnera, mais ici, après de longs préliminaires, le problème est enfin visible - les performances asymptotiquement faibles de ce solveur dans le cas de grands systèmes. Et tout cela à cause du fait d'avoir une boucle for . Quel est le problème? Même si la condition à l'intérieur du if se révèle extrêmement rarement vraie, il n'y a pas d'éloignement du cycle lui-même, ce qui crée la complexité de l'algorithme $$$O (N)$ où $$$N$ Est la taille de la matrice. Peu importe la façon dont les boucles sont optimisées par les compilateurs modernes, cette complexité se fera sentir avec de grandes $$$N$ . Particulièrement offensant lorsque l'ensemble du vecteur $$$f$ se compose entièrement de zéros, car, comme nous l'avons dit, alors ne faites rien dans ce cas! Pas une seule opération arithmétique! Que diable $$$O (N)$ action?

Eh bien, disons. Même dans ce cas, pourquoi ne pouvez-vous pas simplement endurer l'inactivité, car de vrais calculs avec des nombres réels, c'est-à-dire ceux qui relèvent if seront encore peu nombreux? Mais le fait est que cette procédure d'écoulement direct avec un côté droit raréfié lui-même est souvent utilisée dans les cycles externes et sous-tend la décomposition de Cholesky $$$A = LL ^ T$ et décomposition LU à gauche. Oui, oui, une de ces extensions, sans la capacité de faire tous ces mouvements directs et inverses dans la résolution de systèmes linéaires, perd tout sens pratique.

Théorème Si la matrice est définie positive symétrique (SPD), elle peut être représentée comme $$$A = LL ^ T$ la seule façon où $$$L$ - matrice triangulaire inférieure avec des éléments positifs sur la diagonale.

Pour les matrices SPD très clairsemées, une décomposition de Cholesky ascendante est utilisée. Représentation schématique de la décomposition sous forme de bloc matriciel

$$

$$\ begin {pmatrix} \ mathbf {L} _ {11} & \ mathbf {0} \\ \ mathbf {l} _ {12} ^ T & l_ {22} \ end {pmatrix} \ begin {pmatrix} \ mathbf {L} _ {11} ^ T & \ mathbf {l} _ {12} \\ \ mathbf {0} & l_ {22} \ end {pmatrix} = \ begin {pmatrix} \ mathbf {A} _ { 11} & \ mathbf {a} _ {12} \\ \ mathbf {a} _ {12} ^ T & a_ {22} \ end {pmatrix},$$

l'ensemble du processus de factorisation peut être logiquement divisé en seulement trois étapes.

Algorithme 3:

1. méthode supérieure de décomposition de Cholesky de plus petite dimension $$$\ mathbf {L} _ {11} \ mathbf {L} _ {11} ^ T = \ mathbf {A} _ {11}$ (récursivité!)
2. système triangulaire inférieur avec côté droit clairsemé $$$\ mathbf {L} _ {11} \ mathbf {l} _ {12} = \ mathbf {a} _ {12}$ (le voici! )
3. calcul $$$l_ {22} = \ sqrt {a_ {22} - \ mathbf {l} _ {12} ^ T \ mathbf {l} _ {12}}$ (opération triviale)

En pratique, cela est mis en œuvre de telle manière que les étapes 3 et 2 sont effectuées en un seul grand cycle, et dans cet ordre. Ainsi, la matrice est exécutée en diagonale de haut en bas, augmentant la matrice $$$L$ ligne par ligne à chaque itération de la boucle.

Si, dans un tel aglorithme, la complexité de l'étape 2 sera $$$O (N)$ où $$$N$ Est la taille de la matrice triangulaire inférieure $$$\ mathbf {L} _ {11}$ sur une itération arbitraire d'un grand cycle, alors la complexité de la décomposition entière sera au moins $$$O (N ^ 2)$ ! Oh, je ne voudrais pas ça!

Etat de l'art

De nombreux algorithmes sont en quelque sorte basés sur l'imitation d'actions humaines pour résoudre des problèmes. Si vous donnez à une personne un système linéaire triangulaire inférieur avec un côté droit, dans lequel seulement 1-2 sont non nuls, alors bien sûr, il exécute d'abord le vecteur du côté droit avec ses yeux de haut en bas (ce maudit cycle de complexité $$$O (N)$ ! ) pour trouver ces non-zéros. Ensuite, il ne les utilisera que sans perdre de temps sur les composantes nulles, car ces dernières n'affecteront pas la solution: cela n'a aucun sens de diviser les zéros en composantes diagonales de la matrice, ainsi que de déplacer la colonne multipliée par zéro vers la droite. C'est l'algorithme ci-dessus 2. Il n'y a pas de miracles. Mais que se passe-t-il si une personne reçoit immédiatement le nombre de composants non nuls provenant d'autres sources? Par exemple, si le côté droit est une colonne d'une autre matrice, comme c'est le cas avec la décomposition de Cholesky, alors nous avons un accès instantané à ses composants non nuls, si demandé séquentiellement:

 //     j,     : for (size_t p = jA[j]; p < jA[j+1]; p++) //    vA[p]     iA[p]   j 

La complexité de cet accès est $$$O (nnz_j)$ où $$$nnz_j$ Est le nombre de composants différents de zéro dans la colonne j. Dieu merci pour le format CSC! Comme vous pouvez le voir, il n'est pas seulement utilisé pour économiser de la mémoire.

Dans ce cas, nous pouvons saisir l'essence de ce qui se passe lors de la résolution par la méthode du cours direct $$$Lx = f$ pour matrices clairsemées $$$L$ et côté droit $$$f$ . Retenez votre souffle: nous prenons une composante non nulle $$$f_i$ à droite, on trouve la variable correspondante $$$x_i$ diviser par $$$L_ {ii}$ puis, en multipliant la i-ème colonne entière par cette variable trouvée, nous introduisons des non-zéros supplémentaires sur le côté droit en soustrayant cette colonne du côté droit! Ce processus est bien décrit dans le langage des ... graphiques. De plus, orienté et non cyclique.

Pour une matrice triangulaire inférieure qui n'a pas de zéros sur la diagonale, nous définissons un graphe connecté. Nous supposons que la numérotation des lignes et des colonnes commence à zéro.
Définition Un graphique de connectivité pour une matrice triangulaire inférieure $$$L$ la taille $$$N$ , qui n'a pas de zéros sur la diagonale, est appelé l'ensemble des nœuds $$$V = \ {0, ..., N-1 \}$ et bords orientés $$$E = \ {(j, i) | L_ {ij} \ ne 0 \}$ .

Voici, par exemple, à quoi ressemble le graphe de connexion d'une matrice triangulaire inférieure.

dans lequel les nombres correspondent simplement au nombre ordinal de l'élément diagonal, et les points indiquent des composantes non nulles en dessous de la diagonale:



Définition Graphique orienté accessibilité $$$G$ sur plusieurs indices $$$W \ sous-ensemble V$ appeler autant d'indices $$$R_G (W) \ sous-ensemble V$ que dans n'importe quel index $$$z \ dans R_G (W)$ vous pouvez obtenir en passant par le graphique à partir d'un index $$$w (z) \ en W$ .

Exemple: pour un graphique à partir d'une image $$$R_G (\ {0, 4 \}) = \ {0, 4, 5, 6 \}$ . Il est clair que cela tient toujours $$$W \ sous-ensemble R_G (W)$ .

Si nous représentons chaque nœud du graphique comme le numéro de colonne de la matrice qui l'a généré, alors les nœuds voisins vers lesquels ses bords sont dirigés correspondent aux numéros de ligne des composants non nuls de cette colonne.

Soit $$$nnz (x) = \ {j | x_j \ ne 0 \}$ désigne l'ensemble des indices correspondant à des positions non nulles dans le vecteur x.

Théorème de Hilbert (non, pas sous le nom duquel les espaces sont nommés) $$$nnz (x)$ où $$$x$ il y a un vecteur de solution d'un système triangulaire inférieur clairsemé $$$Lx = f$ avec le côté droit clairsemé $$$f$ coïncide avec $$$R_G (nnz (f))$ .

Addition par nous-mêmes: dans le théorème, nous ne prenons pas en compte la possibilité improbable que les nombres non nuls du côté droit, lors de la destruction des colonnes, puissent être réduits à zéro, par exemple, 3 - 3 = 0. Cet effet est appelé annulation numérique. Prendre en compte ces zéros spontanés est une perte de temps et ils sont perçus comme tous les autres nombres dans des positions non nulles.

Une méthode efficace pour effectuer un mouvement direct avec un côté droit clairsemé donné, sous l'hypothèse que nous avons un accès direct à ses composantes non nulles par des indices, est la suivante.

1. Nous parcourons le graphique en utilisant la méthode de la «recherche en profondeur d'abord», en commençant séquentiellement à partir de l'index $$$i \ in nnz (f)$ chaque composant non nul du côté droit. Écriture de nœuds trouvés dans un tableau $$$R_G (nnz (f))$ en faisant cela dans l'ordre dans lequel nous revenons au graphique! Par analogie avec l'armée d'envahisseurs: nous occupons le pays sans cruauté, mais quand nous avons commencé à être repoussés, nous, en colère, détruisons tout sur son passage.
   Il convient de noter qu’il n’est absolument pas nécessaire que la liste des indices $$$nnz (f)$ a été trié par ordre croissant lorsqu'il a été alimenté à l'entrée de l'algorithme de «recherche en profondeur d'abord». Vous pouvez commencer dans n'importe quel ordre sur le plateau $$$nnz (f)$ . Séquence différente appartenant à l'ensemble $$$nnz (f)$ index n'affecte pas la solution finale, comme nous le verrons dans l'exemple.
   
   Cette étape ne nécessite aucune connaissance d'un vrai tableau. $$$v_A$ ! Bienvenue dans le monde de l'analyse symbolique lorsque vous travaillez avec des solveurs directs clairsemés!
   
2. Nous procédons à la solution elle-même, en ayant à notre disposition un tableau $$$R_G (nnz (f))$ de l'étape précédente. Nous détruisons les colonnes dans l' ordre inverse de l'enregistrement de ce tableau . Voila!

Exemple

Prenons un exemple dans lequel les deux étapes sont illustrées en détail. Supposons que nous ayons une matrice 12 x 12 de la forme suivante:

Le graphique de connectivité correspondant a la forme:

Laisser le non nul dans la partie droite uniquement aux positions 3 et 5, c'est-à-dire $$$nnz (f) = \ {3, 5 \}$ . Parcourons le graphique, en commençant par ces indices dans l'ordre écrit. Ensuite, la méthode de "recherche approfondie" se présente comme suit. Nous allons d'abord visiter les indices afin $$ , sans oublier de marquer les index comme visités. Dans la figure ci-dessous, ils sont ombrés en bleu:

Lors du retour, nous entrons des indices dans notre tableau d'indices de composants de solution non nuls $$$nnz (x): = \ {\ color {blue} {11}, \ color {blue} 8, \ color {blue} 3 \}$ . Ensuite, essayez d'exécuter $$$5 \ to8 \ to ...$ , mais nous rencontrons un nœud 8 déjà marqué, nous ne touchons donc pas à cet itinéraire et nous nous rendons à la succursale $$$5 \ to9 \ to ...$ . Le résultat de cette course sera $$ . Nous ne pouvons pas visiter le noeud 11, car il est déjà étiqueté. Alors, revenez en arrière et complétez le tableau $$$nnz (x)$ nouvelle prise marquée en vert: $$$nnz (x): = \ {\ color {blue} {11}, \ color {blue} 8, \ color {blue} 3, \ color {green} {10}, \ color {green} 9, \ color {vert} 5 \}$ . Et voici l'image:

Les nœuds verts boueux 8 et 11 sont ceux que nous voulions visiter lors de la deuxième manche, mais nous n'avons pas pu, car déjà visité lors du premier.

Ainsi, le tableau $$$nnz (x) = R_G (nnz (f))$ formé. Nous passons à la deuxième étape. Une étape simple: nous parcourons le tableau $$$nnz (x)$ dans l'ordre inverse (de droite à gauche), trouver les composants correspondants du vecteur solution $$$x$ division par les composantes diagonales de la matrice $$$L$ et déplacer les colonnes vers la droite. Autres composants $$$x$ comme ils étaient des zéros, ils le sont restés. Schématiquement, cela est illustré ci-dessous, où la flèche du bas indique l'ordre dans lequel les colonnes sont détruites:

Remarque: dans cet ordre de destruction des colonnes, le numéro 3 sera rencontré plus tard que les numéros 5, 9 et 10! Les colonnes ne sont pas détruites dans l'ordre croissant, ce qui serait une erreur pour les colonnes denses, c'est-à-dire matrices inachevées. Mais pour les matrices clairsemées, cela est courant. Comme le montre la structure matricielle non nulle de cet exemple, l'utilisation des colonnes 5, 9 et 10 dans la colonne 3 ne faussera pas les composants dans la réponse $$$x_5$ , $$$x_9$ , $$$x_ {10}$ pas du tout, ils ont c $$$x_3$ juste "pas d'intersections". Notre méthode en a tenu compte! De même, l'utilisation de la colonne 8 après la colonne 9 ne gâchera pas le composant $$$x_9$ . Excellent algorithme, n'est-ce pas?

Si nous parcourons d'abord le graphique à partir de l'index 5, puis à partir de l'index 3, notre tableau sera $$$nnz (x): = \ {\ color {blue} {11}, \ color {blue} 8, \ color {blue} {10}, \ color {blue} {9}, \ color {blue} 5, \ color {green} 3 \}$ . Détruire les colonnes dans l'ordre inverse de ce tableau donnera exactement la même solution que dans le premier cas!

La complexité de toutes nos opérations est mise à l'échelle par le nombre d'opérations arithmétiques réelles et le nombre de composants non nuls sur le côté droit, mais pas par la taille de la matrice! Nous avons atteint notre objectif.

La critique

CEPENDANT! Un lecteur critique remarquera que l'hypothèse même au départ que nous avons «un accès direct aux composants non nuls du côté droit par index» signifie déjà qu'une fois avant, nous avons parcouru le côté droit de haut en bas pour trouver ces mêmes indices et les organiser. dans le tableau $$$nnz (f)$ c'est-à-dire avoir déjà dépensé $$$O (N)$ action! De plus, le graphe exécuté lui-même nécessite que nous allouions d'abord la mémoire de la longueur maximale possible (ailleurs, vous devez écrire les index trouvés en cherchant en profondeur!) Afin de ne pas souffrir d'une nouvelle réaffectation à mesure que le tableau grandit $$$nnz (x)$ , et cela nécessite également $$$O (N)$ opérations! Pourquoi, alors, disent-ils, tout le tapage?

Mais en effet, pour une solution ponctuelle d'un système triangulaire inférieur clairsemé avec un côté droit clairsemé, initialement défini comme un vecteur dense, il n'y a aucun intérêt à perdre du temps de développement sur tous les algorithmes mentionnés ci-dessus. Ils peuvent être encore plus lents que la méthode du front présentée par l'algorithme 2 ci-dessus. Mais, comme déjà mentionné, cet appareil est indispensable pour la factorisation de Cholesky, vous ne devez donc pas me lancer de tomates. En effet, avant d'exécuter l'algorithme 3, toute la mémoire nécessaire de longueur maximale est allouée immédiatement et nécessite $$$O (N)$ le temps. Dans un autre cycle de colonnes $$$A$ toutes les données sont écrasées uniquement dans un tableau de longueur fixe $$$N$ , et uniquement dans les positions où cette réécriture est pertinente, grâce à un accès direct à des éléments non nuls. Et c'est précisément à cause de cela que l'efficacité naît!

Implémentation C ++

Le graphique lui-même en tant que structure de données dans le code n'est pas nécessaire à construire. Il suffit de l'utiliser implicitement lorsque l'on travaille directement avec la matrice. Toute la connectivité requise sera prise en compte par l'algorithme. La recherche en profondeur est bien sûr mise en œuvre de manière pratique en utilisant une récursion banale. Ici, par exemple, ressemble à une recherche de profondeur récursive basée sur un seul index de démarrage:

 /* j -   iA, jA -    ,    CSC top -     nnz(x);     0 result -   N,    nnz(x)   0  top-1  marked -    N;      */ void DepthFirstSearch(size_t j, const std::vector<size_t>& iA, const std::vector<size_t>& jA, size_t& top, std::vector<size_t>& result, std::vector<int>& marked) { size_t p; marked[j] = 1; //   j   for (p = jA[j]; p < jA[j+1]; p++) //       j { if (!marked[iA[p]]) //  iA[p]   { DepthFirstSearch(iA[p], iA, jA, top, result, marked); //      iA[p] } } result[top++] = j; //  j   nnz(x) } 

Si vous passez la variable top égale à zéro dans le tout premier appel à DepthFirstSearch , puis après l'achèvement de la procédure récursive entière, la variable top sera égale au nombre d'index trouvés dans le tableau de result . Devoirs: écrire une autre fonction qui prend un tableau d'indices de composants non nuls sur le côté droit et les transmet séquentiellement à la fonction DepthFirstSearch . Sans cela, l'algorithme n'est pas complet. Veuillez noter: un tableau de nombres réels $$$vA$ nous ne le transmettons pas du tout à la fonction, car il n'est pas nécessaire dans le processus d'analyse symbolique.

Malgré sa simplicité, un défaut d'implémentation est évident: pour les grands systèmes, les débordements de pile sont imminents. Eh bien, il y a une option basée sur une boucle au lieu de la récursivité. C'est plus difficile à comprendre, mais déjà adapté à toutes les occasions:

 /* j -   N -   iA, jA -    ,    CSC top -     nnz(x) result -   N,     nnz(x)   top  ret-1 ,  ret -    marked -    N work_stack -     N */ size_t DepthFirstSearch(size_t j, size_t N, const std::vector<size_t>& iA, const std::vector<size_t>& jA, size_t top, std::vector<size_t>& result, std::vector<int>& marked, std::vector<size_t>& work_stack) { size_t p, head = N - 1; int done; result[N - 1] = j; //    while (head < N) { j = result[head]; //  j     if (!marked[j]) { marked[j] = 1; //   j   work_stack[head] = jA[j]; } done = 1; //    j      for (p = work_stack[head] + 1; p < jA[j+1]; p++) //    j { //   c   iA[p] if (marked[iA[p]]) continue; //  iA[p]  ,   work_stack[head] = p; //       j result[--head] = iA[p]; //       iA[p] done = 0; //   j    break; } if (done) //      j  { head++; //  j    result[top++] = j; //  j    } } return (top); } 

Et voici à quoi ressemble déjà le générateur d'une structure vectorielle de solution non nulle $$$nnz (x)$ :

 /* iA, jA -   ,    CSC iF -        f nnzf -      f result -   N,    nnz(x)   0  ret-1 ,  ret -    marked -    N,    work_stack -     N */ size_t NnzX(const std::vector<size_t>& iA, const std::vector<size_t>& jA, const std::vector<size_t>& iF, size_t nnzf, std::vector<size_t>& result, std::vector<int>& marked, std::vector<size_t>& work_stack) { size_t p, N, top; N = jA.size() - 1; top = 0; for (p = 0; p < nnzf; p++) if (!marked[iF[p]]) //        top = DepthFirstSearch(iF[p], N, iA, jA, vA, top, result, marked, work_stack); for (p = 0; p < top; p++) marked[result[p]] = 0; //   return (top); } 

Enfin, en combinant tout en un, nous écrivons le solveur triangulaire inférieur pour le cas du côté droit raréfié:

 /* iA, jA, vA -  ,    CSC iF -        f nnzf -      f vF -       f result -   N,    nnz(x)   0  ret-1 ,  ret -    marked -    N,    work_stack -     N x -    N,    ;     */ size_t lower_triangular_solve (const std::vector<size_t>& iA, const std::vector<size_t>& jA, const std::vector<double>& vA, const std::vector<size_t>& iF, size_t nnzf, const std::vector<double>& vF, std::vector<size_t>& result, std::vector<int>& marked, std::vector<size_t>& work_stack, std::vector<double>& x) { size_t top, p, j; ptrdiff_t px; top = NnzX(iA, jA, iF, nnzf, result, marked, work_stack); for (p = 0; p < nnzf; p++) x[iF[p]] = vF[p]; //    for (px = top; px > -1; px--) //     { j = result[px]; // x(j)   x [j] /= vA[jA[j]]; //   x(j) for (p = jA[j]+1; p < jA[j+1]; p++) { x[iA[p]] -= vA[p]*x[j]; //  j-  } } return (top) ; } 

On voit que notre cycle ne passe que par les indices du tableau $$$nnz (x)$ , pas l'ensemble $$ . C'est fait!

Il existe une implémentation qui n'utilise pas un tableau de balises markedpour économiser de la mémoire. Au lieu de cela, un ensemble supplémentaire d'indices est utilisé.$$$V_1$ ne se croisant pas avec l'ensemble $$$V = \{0, ..., N - 1\}$dans lequel il y a un mappage un à un par une simple opération algébrique comme procédure de marquage de nœud. Cependant, à notre époque de la mémoire bon marché, l'enregistrer sur un seul tableau de longueur$$$N$ semble complètement redondant.

En conclusion

Le processus de résolution d'un système clairsemé d'équations linéaires par la méthode directe est généralement divisé en trois étapes:

1. Analyse de caractère
2. Factorisation numérique basée sur les données de l'analyse sivol
3. La solution des systèmes triangulaires obtenus avec un côté droit dense

La deuxième étape, la factorisation numérique, est la partie la plus consommatrice de ressources et dévore le plus (> 90%) du temps estimé. Le but de la première étape est de réduire le coût élevé de la seconde. Un exemple d'analyse symbolique a été présenté dans ce billet. Cependant, c'est la première étape qui nécessite le temps de développement le plus long et les coûts mentaux maximaux de la part du développeur. Une bonne analyse symbolique nécessite la connaissance de la théorie des graphes et des arbres et la possession d'un «instinct algorithmique». La deuxième étape est incomparablement plus facile à mettre en œuvre.
Eh bien, la troisième étape, à la fois dans la mise en œuvre et dans le temps estimé, est dans la plupart des cas la plus sans prétention.

Une bonne introduction aux méthodes directes pour les systèmes clairsemés se trouve dans Tim Davis , professeur au Département d'informatique de la Texas A&M University , intitulé "Méthodes directes pour les systèmes linéaires clairsemés ". L'algorithme décrit ci-dessus y est décrit. Je ne connais pas Davis moi-même, même si je travaillais moi-même dans la même université (bien que dans une faculté différente). Si je ne me trompe pas, Davis lui-même a participé au développement de solveurs (!). par ailleurs utilisé dans Matlab, il a été impliqué , même pour les producteurs d'électricité image dans Google Street View (OLS) par ailleurs, à la même faculté énumérés nul autre que lui - même Bjarne Stroustrup - .. le créateur de C ++

peut-être un jour Je vais écrire autre chose sur les matrices clairsemées ou les méthodes numériques dans général. Bon à tous!