Vorspiel

Die numerische Lösung linearer Gleichungssysteme ist in vielen Bereichen der angewandten Mathematik, des Ingenieurwesens und der IT-Branche ein unverzichtbarer Schritt, unabhängig davon, ob sie mit Grafiken arbeitet, die Aerodynamik eines Flugzeugs berechnet oder die Logistik optimiert. Die jetzt modische „Maschine“ wäre auch ohne sie nicht viel vorangekommen. Darüber hinaus macht die Lösung linearer Systeme in der Regel den größten Prozentsatz aller Rechenkosten aus. Es ist klar, dass diese kritische Komponente eine Optimierung der maximalen Geschwindigkeit erfordert.

Arbeiten Sie oft mit den sogenannten spärliche Matrizen - solche, deren Nullen um Größenordnungen größer sind als andere Werte. Dies ist beispielsweise unvermeidlich, wenn Sie sich mit partiellen Differentialgleichungen oder anderen Prozessen befassen, bei denen die entstehenden Elemente in den definierenden linearen Beziehungen nur mit „Nachbarn“ verbunden sind. Hier ist ein mögliches Beispiel einer spärlichen Matrix für die in der klassischen Physik bekannte eindimensionale Poisson-Gleichung $$$- \ phi ^ {''} = f$ auf einem einheitlichen Gitter (ja, es sind zwar nicht so viele Nullen darin, aber beim Schleifen des Gitters sind sie gesund!):

$$

$$\ 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}$$

Die ihnen gegenüberliegenden Matrizen - diejenigen, bei denen sie die Anzahl der Nullen nicht beachten und ausnahmslos alle Komponenten berücksichtigen - werden als dicht bezeichnet.

Sparse-Matrizen werden häufig aus verschiedenen Gründen in einem komprimierten Spaltenformat dargestellt - Compressed Sparse Column (CSC). Diese Beschreibung verwendet zwei Integer-Arrays und ein Gleitkomma. Lassen Sie die Matrix alles haben $$$nnz_A$ ungleich Null und $$$N$ Spalten. Matrixelemente werden ohne Ausnahmen in Spalten von links nach rechts aufgelistet. Erstes Array $$$i_A$ Längen $$$nnz_A$ enthält Zeilennummern von Matrixkomponenten ungleich Null. Zweites Array $$$j_A$ Längen $$$N + 1$ enthält ... nein, nicht die Spaltennummern, da dann die Matrix im Koordinatenformat (Koordinate) oder ternär (Triplett) geschrieben würde. Das zweite Array enthält die Seriennummern der Matrixkomponenten, mit denen die Spalten beginnen, einschließlich einer zusätzlichen Dummy-Spalte am Ende. Schließlich das dritte Array $$$v_A$ Längen $$$nnz_A$ enthält bereits die Komponenten selbst in der Reihenfolge, die dem ersten Array entspricht. Hier beispielsweise unter der Annahme, dass die Nummerierung von Zeilen und Spalten von Matrizen für eine bestimmte Matrix bei Null beginnt

$$

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

Arrays werden sein $$$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 \}$ .

Methoden zur Lösung linearer Systeme werden in zwei große Klassen unterteilt - direkt und iterativ. Die geraden Linien basieren auf der Möglichkeit, die Matrix als Produkt zweier einfacherer Matrizen darzustellen, um die Lösung dann in zwei einfache Stufen aufzuteilen. Iterativ nutzen Sie die einzigartigen Eigenschaften linearer Räume und arbeiten an der Antike als Weltmethode für die sukzessive Approximation von Unbekannten an die gewünschte Lösung. Im Konvergenzprozess selbst werden die Matrizen in der Regel nur verwendet, um sie mit Vektoren zu multiplizieren. Iterative Methoden sind viel billiger als direkte, funktionieren jedoch bei schlecht konditionierten Matrizen nur schleppend. Wenn die Zuverlässigkeit von Stahlbeton wichtig ist, verwenden Sie gerade Linien. Hier bin ich hier und möchte ein wenig ansprechen.

Sagen wir für eine quadratische Matrix $$$A$ Die Zerlegung der Form steht uns zur Verfügung $$$A = LU$ wo $$$L$ und $$$U$ - jeweils die untere Dreiecks- und die obere Dreiecksmatrize. Das erste bedeutet, dass es eine Null über der Diagonale hat, das zweite bedeutet, dass es unter der Diagonale liegt. Wie genau wir diese Zersetzung bekommen haben - wir sind jetzt nicht interessiert. Hier ist ein einfaches Zerlegungsbeispiel:

$$

$$\ 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}$$

Wie in diesem Fall das Gleichungssystem zu lösen $$$Ax = f$ Zum Beispiel mit der rechten Seite $$$f = \ begin {pmatrix} 4 \\ 2 \\ 3 \ end {pmatrix}$ ? Die erste Stufe ist eine direkte Bewegung (Vorwärtslösung = Vorwärtssubstitution). Wir bezeichnen $$$y: = Ux$ und mit dem System arbeiten $$$Ly = f$ . Weil $$$L$ - unteres Dreieck, sukzessive im Zyklus finden wir alle Komponenten $$$y$ von oben nach unten:

$$

$$\ 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} \ impliziert y = \ begin {pmatrix} 4 \\ -6 \\ -1 \ end {pmatrix}$$

Die zentrale Idee ist, dass nach dem Finden $$$i$ Vektorkomponenten $$$y$ Es wird mit einer Spalte mit derselben Matrixnummer multipliziert $$$L$ welches dann von der rechten Seite abgezogen wird. Die Matrix selbst scheint von links nach rechts zu kollabieren und nimmt ab, wenn immer mehr Vektorkomponenten gefunden werden $$$y$ . Dieser Vorgang wird als Spalteneliminierung bezeichnet.

Die zweite Stufe ist eine Rückwärtsbewegung (Rückwärtslösung = Rückwärtssubstitution). Einen Vektor gefunden haben $$$y$ wir entscheiden $$$Ux = y$ . Hier gehen wir die Komponenten bereits von unten nach oben durch, aber die Idee bleibt dieselbe: $$$i$ Die Spalte wird mit der gerade gefundenen Komponente multipliziert $$$x_i$ und bewegt sich nach rechts, und die Matrix kollabiert von rechts nach links. Der gesamte Prozess ist für die im Beispiel im folgenden Bild erwähnte Matrix dargestellt:

$$

$$\ 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} \ impliziert \ begin {pmatrix} 1 & 0 \\ 2 & 1 \ end {pmatrix} \ begin {pmatrix} y_2 \\ y_3 \ end {pmatrix } = \ begin {pmatrix} -6 \\ -13 \ end {pmatrix} \ impliziert \ 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} \ impliziert \ begin {pmatrix} 1 \ end {pmatrix} \ begin {pmatrix} x_1 \ end {pmatrix} = \ begin {pmatrix} -3 \ end {pmatrix}$$

und unsere Entscheidung wird sein $$$x = \ begin {pmatrix} -3 \\ - 9 \\ 2 \ end {pmatrix}$ .

Wenn die Matrix dicht ist, dh vollständig in Form eines eindimensionalen oder zweidimensionalen Arrays dargestellt wird, erfolgt der Zugriff auf ein bestimmtes Element im Laufe der Zeit $$$O (1)$ Dann ist ein ähnliches Lösungsverfahren mit der vorhandenen Zerlegung nicht schwierig und leicht zu codieren, sodass wir nicht einmal Zeit damit verschwenden. Was ist, wenn die Matrix dünn ist? Auch hier ist das grundsätzlich nicht schwierig. Hier ist der C ++ - Code für die Vorwärtsbewegung, in der sich die Lösung befindet $$$x$ Es wird für die Stelle auf der rechten Seite geschrieben, ohne die Richtigkeit der Eingabedaten zu überprüfen (CSC-Arrays entsprechen der unteren Dreiecksmatrix):

Algorithmus 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] ; } } 

Alle weiteren Diskussionen betreffen nur den Vorwärtshub zum Lösen des unteren Dreieckssystems $$$Lx = f$ .

Krawatte

Aber was ist, wenn die rechte Seite, d.h. Vektor rechts vom Gleichheitszeichen $$$Lx = f$ Hat es selbst eine große Anzahl von Nullen? Dann ist es sinnvoll, die mit Nullpositionen verbundenen Berechnungen zu überspringen. Änderungen im Code sind in diesem Fall minimal:

Algorithmus 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]; } } } 

Das einzige, was wir hinzugefügt haben, ist die if , deren Zweck darin besteht, die Anzahl der arithmetischen Operationen auf ihre tatsächliche Anzahl zu reduzieren. Wenn zum Beispiel die gesamte rechte Seite aus Nullen besteht, müssen Sie nichts zählen: Die Lösung ist die rechte Seite.

Alles sieht gut aus und natürlich wird es funktionieren, aber hier ist nach einem langen Vorspiel endlich das Problem sichtbar - die asymptotisch niedrige Leistung dieses Lösers bei großen Systemen. Und das alles wegen der Tatsache, dass es eine for Schleife gibt. Was ist das Problem? Selbst wenn sich herausstellt, dass die Bedingung innerhalb des if äußerst selten wahr ist, gibt es kein Entkommen aus dem Zyklus selbst, und dies schafft die Komplexität des Algorithmus $$$O (N)$ wo $$$N$ Ist die Größe der Matrix. Unabhängig davon, wie die Loops von modernen Compilern optimiert werden, wird sich diese Komplexität mit großem Aufwand bemerkbar machen $$$N$ . Besonders anstößig, wenn der ganze Vektor $$$f$ besteht ganz aus Nullen, denn dann, wie gesagt, dann in diesem Fall nichts tun! Keine einzige arithmetische Operation! Was zur Hölle $$$O (N)$ Aktion?

Sagen wir mal. Warum können Sie den Lauf im Leerlauf nicht einfach ertragen, weil echte Berechnungen mit reellen Zahlen, d. H. diejenigen, die unter fallen, if noch wenige sein werden? Tatsache ist jedoch, dass dieses Direktflussverfahren mit einer verdünnten rechten Seite selbst häufig in externen Zyklen angewendet wird und der Cholesky-Zersetzung zugrunde liegt $$$A = LL ^ T$ und links aussehende LU-Zerlegung . Ja, ja, eine dieser Erweiterungen, ohne die Fähigkeit, all diese direkten und umgekehrten Bewegungen beim Lösen linearer Systeme auszuführen, verliert jede praktische Bedeutung.

Satz Wenn die Matrix symmetrisch positiv definit (SPD) ist, kann sie als dargestellt werden $$$A = LL ^ T$ der einzige Weg wo $$$L$ - untere dreieckige Matrix mit positiven Elementen auf der Diagonale.

Für sehr spärliche SPD-Matrizen wird eine nach oben gerichtete Cholesky-Zersetzung verwendet. Schematische Darstellung der Zerlegung in Matrixblockform

$$

$$\ 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},$$

Der gesamte Faktorisierungsprozess kann logisch in nur drei Schritte unterteilt werden.

Algorithmus 3:

1. obere Methode der Cholesky-Zersetzung kleinerer Dimension $$$\ mathbf {L} _ {11} \ mathbf {L} _ {11} ^ T = \ mathbf {A} _ {11}$ (Rekursion!)
2. unteres Dreieckssystem mit spärlicher rechter Seite $$$\ mathbf {L} _ {11} \ mathbf {l} _ {12} = \ mathbf {a} _ {12}$ ( hier ist es! )
3. Berechnung $$$l_ {22} = \ sqrt {a_ {22} - \ mathbf {l} _ {12} ^ T \ mathbf {l} _ {12}}$ (triviale Operation)

In der Praxis wird dies so implementiert, dass die Schritte 3 und 2 in einem großen Zyklus und in dieser Reihenfolge ausgeführt werden. Somit wird die Matrix diagonal von oben nach unten verlaufen, wodurch die Matrix vergrößert wird $$$L$ Zeile für Zeile bei jeder Iteration der Schleife.

Wenn in einem solchen Aglorithmus die Vorwärtskomplexität in Schritt 2 sein wird $$$O (N)$ wo $$$N$ Ist die Größe der unteren Dreiecksmatrix $$$\ mathbf {L} _ {11}$ Bei einer beliebigen Iteration eines großen Zyklus beträgt die Komplexität der gesamten Zerlegung mindestens $$$O (N ^ 2)$ ! Oh, das würde mir nicht gefallen!

Stand der Technik

Viele Algorithmen basieren irgendwie darauf, menschliche Handlungen bei der Lösung von Problemen nachzuahmen. Wenn Sie einer Person ein unteres dreieckiges lineares System mit einer rechten Seite geben, in der nur 1-2 ungleich Null sind, dann lässt sie natürlich zuerst den Vektor der rechten Seite mit den Augen von oben nach unten laufen (dieser verdammte Zyklus der Komplexität) $$$O (N)$ ! ), um diese Nichtzeros zu finden. Dann wird er nur sie verwenden, ohne Zeit für die Nullkomponenten zu verschwenden, da letztere die Lösung nicht beeinflussen: Es ist nicht sinnvoll, Nullen in die diagonalen Komponenten der Matrix zu unterteilen und die mit Null multiplizierte Spalte nach rechts zu verschieben. Dies ist der obige Algorithmus 2. Es gibt keine Wunder. Aber was ist, wenn eine Person sofort die Anzahl der Komponenten ungleich Null aus anderen Quellen erhält? Wenn zum Beispiel die rechte Seite eine Spalte einer anderen Matrix ist, wie dies bei Choleskys Zerlegung der Fall ist, haben wir dort sofort Zugriff auf die Komponenten ungleich Null, wenn dies nacheinander angefordert wird:

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

Die Komplexität dieses Zugriffs ist $$$O (nnz_j)$ wo $$$nnz_j$ Ist die Anzahl der Nicht-Null-Komponenten in Spalte j. Gott sei Dank für das CSC-Format! Wie Sie sehen, wird es nicht nur zum Speichern von Speicher verwendet.

In diesem Fall können wir die Essenz dessen erfassen, was passiert, wenn wir mit der direkten Kursmethode lösen $$$Lx = f$ für spärliche Matrizen $$$L$ und rechte Seite $$$f$ . Halten Sie den Atem an: Wir nehmen eine Nicht-Null-Komponente $$$f_i$ Auf der rechten Seite finden wir die entsprechende Variable $$$x_i$ Teilen durch $$$L_ {ii}$ und dann, indem wir die gesamte i-te Spalte mit dieser gefundenen Variablen multiplizieren, führen wir zusätzliche Nicht-Nullen auf der rechten Seite ein, indem wir diese Spalte von der rechten Seite subtrahieren! Dieser Prozess ist in der Sprache der ... Grafiken gut beschrieben. Darüber hinaus orientiert und nicht zyklisch.

Für eine untere Dreiecksmatrix ohne Nullen in der Diagonale definieren wir einen verbundenen Graphen. Wir gehen davon aus, dass die Nummerierung von Zeilen und Spalten bei Null beginnt.
Definition Ein Konnektivitätsdiagramm für eine untere Dreiecksmatrix $$$L$ die Größe $$$N$ wird als Knotenmenge bezeichnet, die in der Diagonale keine Nullen aufweist $$$V = \ {0, ..., N-1 \}$ und orientierte Kanten $$$E = \ {(j, i) | L_ {ij} \ ne 0 \}$ .

So sieht beispielsweise der Verbindungsgraph für eine untere Dreiecksmatrix aus.

wobei die Zahlen einfach der Ordnungszahl des diagonalen Elements entsprechen und die Punkte Nicht-Null-Komponenten unterhalb der Diagonale angeben:



Definition Erreichbarkeitsorientiertes Diagramm $$$G$ auf mehreren Indizes $$$W \ Teilmenge V$ nenne so viele Indizes $$$R_G (W) \ Teilmenge V$ das in jedem Index $$$z \ in R_G (W)$ Sie können durch Durchlaufen des Diagramms von einem Index erhalten $$$w (z) \ in W$ .

Beispiel: für eine Grafik aus einem Bild $$$R_G (\ {0, 4 \}) = \ {0, 4, 5, 6 \}$ . Es ist klar, dass es immer gilt $$$W \ Teilmenge R_G (W)$ .

Wenn wir jeden Knoten des Diagramms als die Spaltennummer der Matrix darstellen, die ihn generiert hat, entsprechen die benachbarten Knoten, in die seine Kanten gerichtet sind, den Zeilennummern der Komponenten ungleich Null in dieser Spalte.

Lass $$$nnz (x) = \ {j | x_j \ ne 0 \}$ bezeichnet den Satz von Indizes, die Positionen ungleich Null im Vektor x entsprechen.

Hilberts Satz (nein, nicht nach dessen Namen die Leerzeichen benannt sind) $$$nnz (x)$ wo $$$x$ Es gibt einen Lösungsvektor eines spärlichen unteren Dreieckssystems $$$Lx = f$ mit spärlicher rechter Seite $$$f$ fällt mit zusammen $$$R_G (nnz (f))$ .

Addition allein: Im Theorem berücksichtigen wir nicht die unwahrscheinliche Möglichkeit, dass die Zahlen ungleich Null auf der rechten Seite beim Zerstören der Spalten auf Null reduziert werden können, z. B. 3 - 3 = 0. Dieser Effekt wird als numerische Löschung bezeichnet. Solche spontan auftretenden Nullen zu berücksichtigen, ist Zeitverschwendung und wird als alle anderen Zahlen in Positionen ungleich Null wahrgenommen.

Eine effektive Methode zur Durchführung einer direkten Bewegung mit einer bestimmten spärlichen rechten Seite unter der Annahme, dass wir über Indizes direkten Zugriff auf ihre Nicht-Null-Komponenten haben, ist wie folgt.

1. Wir durchlaufen das Diagramm mit der Methode "Tiefensuche", beginnend nach dem Index $$$i \ in nnz (f)$ jede Nicht-Null-Komponente der rechten Seite. Gefundene Knoten in ein Array schreiben $$$R_G (nnz (f))$ Dabei in der Reihenfolge, in der wir zum Diagramm zurückkehren! In Analogie zur Armee der Invasoren: Wir besetzen das Land ohne Grausamkeit, aber als wir zurückgedrängt wurden, zerstören wir wütend alles auf seinem Weg.
   Es ist erwähnenswert, dass es absolut nicht notwendig ist, die Liste der Indizes $$$nnz (f)$ wurde nach aufsteigend sortiert, als es der Eingabe des Algorithmus "Tiefensuche" zugeführt wurde. Sie können am Set in beliebiger Reihenfolge beginnen $$$nnz (f)$ . Unterschiedliche Reihenfolge zum Set $$$nnz (f)$ Indizes wirken sich nicht auf die endgültige Lösung aus, wie wir im Beispiel sehen werden.
   
   Dieser Schritt erfordert keine Kenntnis eines realen Arrays. $$$v_A$ ! Willkommen in der Welt der symbolischen Analyse bei der Arbeit mit direkten, spärlichen Lösern!
   
2. Wir fahren mit der Lösung selbst fort und verfügen über ein Array $$$R_G (nnz (f))$ aus dem vorherigen Schritt. Wir zerstören die Spalten in umgekehrter Reihenfolge des Datensatzes dieses Arrays . Voila!

Beispiel

Stellen Sie sich ein Beispiel vor, in dem beide Schritte ausführlich demonstriert werden. Angenommen, wir haben eine 12 x 12-Matrix der folgenden Form:

Das entsprechende Konnektivitätsdiagramm hat die Form:

Der Wert ungleich Null im rechten Teil sei nur an den Positionen 3 und 5, d.h. $$$nnz (f) = \ {3, 5 \}$ . Lassen Sie uns das Diagramm ausgehend von diesen Indizes in der schriftlichen Reihenfolge durchgehen. Dann sieht die Methode der "tiefen Suche" wie folgt aus. Zuerst werden wir die Indizes der Reihe nach besuchen $$$3 \ bis 8 \ bis 11$ Nicht zu vergessen, Indizes als besucht zu markieren. In der folgenden Abbildung sind sie blau schattiert:

Bei der Rückkehr geben wir Indizes in unser Array von Indizes von Lösungskomponenten ungleich Null ein $$$nnz (x): = \ {\ color {blue} {11}, \ color {blue} 8, \ color {blue} 3 \}$ . Versuchen Sie als nächstes zu laufen $$$5 \ to8 \ to ...$ Wir stoßen jedoch auf einen bereits markierten Knoten 8, sodass wir diese Route nicht berühren und zum Zweig gehen $$$5 \ to9 \ to ...$ . Das Ergebnis dieses Laufs wird sein $$$5 \ to9 \ to10$ . Wir können Knoten 11 nicht besuchen, da er bereits beschriftet ist. Gehen Sie also zurück und ergänzen Sie das Array $$$nnz (x)$ neuer Fang grün markiert: $$$nnz (x): = \ {\ color {blue} {11}, \ color {blue} 8, \ color {blue} 3, \ color {green} {10}, \ color {green} 9, \ color {grün} 5 \}$ . Und hier ist das Bild:

Die schlammigen grünen Knoten 8 und 11 wollten wir im zweiten Lauf besuchen, konnten es aber nicht, weil schon während des ersten besucht.

Also das Array $$$nnz (x) = R_G (nnz (f))$ gebildet. Wir gehen zum zweiten Schritt über. Ein einfacher Schritt: Wir durchlaufen das Array $$$nnz (x)$ in umgekehrter Reihenfolge (von rechts nach links) die entsprechenden Komponenten des Lösungsvektors finden $$$x$ Division durch die diagonalen Komponenten der Matrix $$$L$ und Verschieben der Spalten nach rechts. Andere Komponenten $$$x$ da sie Nullen waren, blieben sie es auch. Dies ist schematisch unten dargestellt, wobei der untere Pfeil die Reihenfolge angibt, in der die Spalten zerstört werden:

Hinweis: In dieser Reihenfolge der Zerstörung von Spalten wird Nummer 3 später als Nummer 5, 9 und 10 angetroffen! Spalten werden nicht in aufsteigender Reihenfolge zerstört, was ein Fehler für dichte wäre, d. H. unvollständige Matrizen. Bei spärlichen Matrizen ist dies jedoch üblich. Wie aus der Nicht-Null-Matrixstruktur in diesem Beispiel ersichtlich ist, verzerrt die Verwendung der Spalten 5, 9 und 10 bis Spalte 3 die Komponenten in der Antwort nicht $$$x_5$ , $$$x_9$ , $$$x_ {10}$ überhaupt nicht, sie haben c $$$x_3$ nur "keine Kreuzungen." Unsere Methode hat dies berücksichtigt! In ähnlicher Weise wird die Komponente durch die Verwendung von Spalte 8 nach Spalte 9 nicht beeinträchtigt $$$x_9$ . Großartiger Algorithmus, nicht wahr?

Wenn wir zuerst von Index 5 und dann von Index 3 um den Graphen herumgehen, ist unser Array $$$nnz (x): = \ {\ color {blue} {11}, \ color {blue} 8, \ color {blue} {10}, \ color {blue} {9}, \ color {blue} 5, \ color {green} 3 \}$ . Wenn Sie die Spalten in umgekehrter Reihenfolge dieses Arrays zerstören, erhalten Sie genau die gleiche Lösung wie im ersten Fall!

Die Komplexität aller unserer Operationen wird durch die Anzahl der tatsächlichen arithmetischen Operationen und die Anzahl der Komponenten ungleich Null auf der rechten Seite skaliert, jedoch nicht durch die Größe der Matrix! Wir haben unser Ziel erreicht.

Kritik

JEDOCH! Ein kritischer Leser wird feststellen, dass die Annahme zu Beginn, dass wir „direkten Zugriff auf die Nicht-Null-Komponenten der rechten Seite über den Index“ haben, bereits einmal bedeutet, dass wir die rechte Seite einmal von oben nach unten durchlaufen haben, um genau diese Indizes zu finden und zu organisieren im Array $$$nnz (f)$ , das heißt, haben bereits ausgegeben $$$O (N)$ Aktion! Darüber hinaus erfordert der Graph-Lauf selbst, dass wir zuerst den Speicher mit der maximal möglichen Länge zuweisen (an einer anderen Stelle müssen Sie die durch eingehende Suche gefundenen Indizes notieren!), Um nicht an einer weiteren Neuzuweisung zu leiden, wenn das Array wächst $$$nnz (x)$ , und das erfordert auch $$$O (N)$ Operationen! Warum dann, sagen sie, die ganze Aufregung?

In der Tat ist es für eine einmalige Lösung eines spärlichen unteren Dreieckssystems mit einer spärlichen rechten Seite, die ursprünglich als dichter Vektor definiert wurde, nicht sinnvoll, Entwicklerzeit für alle oben genannten Algorithmen zu verschwenden. Sie können sogar langsamer sein als die Stirnmethode, die durch den obigen Algorithmus 2 dargestellt wird. Aber wie bereits erwähnt, ist dieser Apparat für die Faktorisierung von Cholesky unverzichtbar, deshalb sollten Sie keine Tomaten auf mich werfen. In der Tat wird vor dem Ausführen von Algorithmus 3 der gesamte erforderliche Speicher mit maximaler Länge sofort zugewiesen und benötigt $$$O (N)$ Zeit. In einem weiteren Spaltenzyklus $$$A$ Alle Daten werden nur in einem Array fester Länge überschrieben $$$N$ und nur an den Positionen, an denen diese Umschreibung relevant ist, dank des direkten Zugriffs auf Elemente ungleich Null. Und genau deshalb entsteht Effizienz!

C ++ Implementierung

Das Diagramm selbst als Datenstruktur im Code muss nicht erstellt werden. Es reicht aus, es implizit zu verwenden, wenn Sie direkt mit der Matrix arbeiten. Alle erforderlichen Konnektivitäten werden vom Algorithmus berücksichtigt. Die Tiefensuche wird natürlich bequem mit banaler Rekursion durchgeführt. Hier sieht es beispielsweise wie eine rekursive Tiefensuche aus, die auf einem einzelnen Startindex basiert:

 /* 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) } 

Wenn Sie beim ersten Aufruf von DepthFirstSearch die top Variable gleich Null übergeben, entspricht die top Variable nach Abschluss der gesamten rekursiven Prozedur der Anzahl der im result gefundenen Indizes. Hausaufgabe: Schreiben Sie eine weitere Funktion, die ein Array von Indizes von Komponenten ungleich Null auf der rechten Seite verwendet und diese nacheinander an die DepthFirstSearch Funktion DepthFirstSearch . Ohne dies ist der Algorithmus nicht vollständig. Bitte beachten Sie: eine Reihe von reellen Zahlen $$$vA$ Wir geben es überhaupt nicht an die Funktion weiter, weil es wird im Prozess der symbolischen Analyse nicht benötigt.

Trotz seiner Einfachheit ist ein Fehler in der Implementierung offensichtlich: Bei großen Systemen stehen Stapelüberläufe vor der Tür. Dann gibt es eine Option, die auf einer Schleife anstelle einer Rekursion basiert. Es ist schwieriger zu verstehen, aber bereits für alle Gelegenheiten geeignet:

 /* 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); } 

Und so sieht der Generator einer Vektorstruktur ungleich Null bereits aus $$$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); } 

Schließlich schreiben wir, indem wir alles zu einem kombinieren, den unteren Dreieckslöser für den Fall der verdünnten rechten Seite:

 /* 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) ; } 

Wir sehen, dass unser Zyklus nur durch die Indizes des Arrays läuft $$$nnz (x)$ , nicht das ganze Set $$$0, 1, ..., N - 1$ . Fertig!

Es gibt eine Implementierung, die kein Array von Tags markedverwendet, um Speicherplatz zu sparen. Stattdessen wird ein zusätzlicher Satz von Indizes verwendet.$$$V_1$ nicht mit der Menge schneiden $$$V = \{0, ..., N - 1\}$in die es eine Eins-zu-Eins-Zuordnung durch eine einfache algebraische Operation als Knotenmarkierungsprozedur gibt. In unserer Zeit des billigen Speichers sparen Sie ihn jedoch auf einem einzigen Array von Längen$$$N$ scheint völlig überflüssig.

Abschließend

Der Prozess der Lösung eines spärlichen linearen Gleichungssystems nach der direkten Methode ist in der Regel in drei Stufen unterteilt:

1. Charakteranalyse
2. Numerische Faktorisierung basierend auf den Daten der Sivol-Analyse
3. Die Lösung der erhaltenen Dreieckssysteme mit einer dichten rechten Seite

Der zweite Schritt, die numerische Faktorisierung, ist der ressourcenintensivste Teil und verschlingt den größten Teil (> 90%) der geschätzten Zeit. Der Zweck des ersten Schritts besteht darin, die hohen Kosten des zweiten zu reduzieren. In diesem Beitrag wurde ein Beispiel für eine symbolische Analyse vorgestellt. Dies ist jedoch der erste Schritt, der die längste Entwicklungszeit und maximale mentale Kosten seitens des Entwicklers erfordert. Eine gute symbolische Analyse erfordert Kenntnisse der Theorie von Graphen und Bäumen und den Besitz eines „algorithmischen Instinkts“. Der zweite Schritt ist unvergleichlich einfacher zu implementieren.
Nun, der dritte Schritt, sowohl bei der Implementierung als auch in der geschätzten Zeit, ist in den meisten Fällen der unprätentiöseste.

Eine gute Einführung in direkte Methoden für spärliche Systeme bietet Tim Davis , Professor am Institut für Informatik der Texas A & M University mit dem Titel "Direkte Methoden für spärliche lineare Systeme ". Der oben beschriebene Algorithmus ist dort beschrieben. Ich kenne Davis selbst nicht, obwohl ich selbst an derselben Universität gearbeitet habe (wenn auch an einer anderen Fakultät). Wenn ich mich nicht irre, war Davis selbst an der Entwicklung von Lösern beteiligt (!). verwendet darüber hinaus in Matlab, er wurde sogar für Stromgeneratoren Bild in Google Street View beteiligt (OLS) durch die Art und Weise, in der gleichen Fakultät aufgeführt nichts anderes als sich selbst Bjarne Stroustrup - .. der Schöpfer von C ++

Vielleicht eines Tages Ich werde noch etwas über spärliche Matrizen oder numerische Methoden schreiben allgemein. Gut zu allen!