nicht gleich an, dass man sie durch Umordnen von Zeilen und Spalten in die tridiagonale Matrix (das ist eine (1,1)-Bandmatrix)
A˜:=
4 1 0 0 0 0 0 1 4 1 0 0 0 0 0 1 4 1 0 0 0 0 0 1 4 1 0 0 0 0 0 1 4 1 0 0 0 0 0 1 4 1 0 0 0 0 0 1 4
umwandeln kann.
Wir wollen uns deshalb im nächsten Unterabschnitt mit dem Problem auseinandersetzen, die Zeilen und Spalten einer gegebenen Matrix so zu vertauschen, dass die Matrix eine Bandmatrix möglichst kleiner Breite wird.
Anmerkungen 4.5
1. Die Ecke vj des Graphen ist sowohl der j-ten Variable als auch der j-ten Gleichung zugeordnet. Der Matrix
B =
5 2 0 0 0 −1
2 4 0 3 0 0
0 0 3 0 0 −2
0 3 0 4 0 3
0 0 0 0 7 1
−1 0 −2 3 1 20
.
ist so z.B. der Graph
0 0.5 1 1.5 2 2.5 3 3.5 4
0 0.5 1 1.5 2 2.5 3 3.5 4
Adjazenzgraph zu B
1 2
3 6 4
5
Abbildung 82: Adjazenzgraph zur Matrix B
zugeordnet. Die blauen Zahlen an den Ecken sind die Nummern j der vj. Knoten k ist mit Knoten i verbunden, wenn aij = aji ̸= 0. Eine Umnumerierung der Ecken entspricht einer simultanen Vertauschung von Variablen und zugehörigen Gleichun-gen:
A−→P APT, P = Permutationsmatrix
2. Die Gleichungsstruktur nichtsymmetrischer Gleichungssysteme modelliert man z.B.
durch gerichtete Graphen. Dabei gehört das geordnete Paar(vi, vj)zur Kantenmenge von G(A), wennaij ̸= 0.(vj, vi)gehört nur dann auch dazu, wenn auch aj,i ̸= 0. Wir gehen auf den nichtsymmetrischen Fall hier nicht ein.
3. In der Graphentheorie verbindet man umgekehrt mit einem ungerichteten Graphen eine symmetrische Adjazenzmatrix A, die gerade für jede Kante {vi, vj} die Ele-mente aij und aji gleich 1 setzt und den anderen den Wert Null gibt. Für gerichtete Graphen wählt man analog eine (möglicherweise) unsymmetrischeAdjazenzmatrix.
4. Diskretisierungsstrukturen tragen oft selbst schon die Struktur des Adjazenz-Graphen.
Das Gitter aus Abbildung 7 ist so z.B. der Adjazenzgraph der auf das Bild folgenden Matrix A. Ersetzt man umgekehrt in A alle von Null verschiedenen Werte durch 1, erhält man die Adjazenzmatrix zum Gittergraphen.
Das Ziel ist nun, eine solche Umnumerierung der Ecken eines Adjazenzgraphen zu finden, bei denen im Netz benachbarte Ecken möglichst benachbarte Eckennummern erhalten.
Leider ist das Problem, die optimale Numerierung zu finden171, NP-schwer172. Daher ist man bei seiner Behandlung auf heuristische Algorithmen angewiesen.
Ein einfaches — aber erfolgreiches — Exemplar solcher heuristischer Algorithmen ist der sogenannte Algorithmus von Cuthill und McKee.
Zu seiner Formulierung benötigen wir nur noch einige Sprechweisen für die Ecken und Kanten des Adjazenzgraphen einer Matrix:
Sprechweisen:
a. vi ist Nachbar vonvj, wenn{vi, vj} ∈K.(auf hochdeutsch: wenn sie verbunden sind) b. Eine Kante k inzidiert mit einer Ecke vi, wenn vi ∈ k ist. (hochdeutsch: wenn die
Ecke von der Kante getroffen wird)
c. Der Grad einer Eckevi ist die Anzahl der mit ihr inzidierenden verschiedenen Kanten.
( = Anzahl der Kanten, dievi mit anderen Ecken verbinden)
Damit können wir den Cuthill-McKee- oder kürzer den CM-Algorithmus nur einfach for-mulieren. In der Graphentheorie ist er auch ùnter dem Namen „Breitensuche“ bekannt173.
Algorithmus 4.6 (Algorithmus von Cuthill-McKee.) Finde Startecke v1 (Wie? Siehe unten.)
For i:= 1 to n−1
Finde alle unnumerierten Nachbarn vonvi und numeriere sie aufsteigend nach steigendem Grad,
beginnend mit der nächsten noch nicht vergebenen Nummer.
end
Bemerkung: Damit der Algorithmus den ganzen Graphen durchnumerieren kann, müs-sen wir annehmen, dass der Adjazenzgraph des behandelten Systems zusammenhängend ist (d.h. dass er nicht in zwei oder mehr nicht durch Kanten verbundene Teilgraphen zer-fällt). Andernfalls bräche der Algorithmus nach Durchnumerierung des Teilgraphen, derv1 enthält (der sogenannten Zusammenhangskomponente von v1), ab.
Das wäre aber nicht schlimm. Man müsste den Algorithmus einfach nur mit einer noch nicht numerierten Ecke neu starten und erhielte am Ende alle Zusammenhangskompo-nenten in durchnumerierter Form. Tatsächlich hätte man damit das Gleichungssytem in ebenso viele unabhängige Gleichungssysteme separiert, die alle einzeln weiterbehandelt werden könnten.
Beispiel CM1: Wir wollen den CM-Algorithmus an einem Beispiel durchspielen. Ge-geben sei dazu die folgende Diskretisierungsstruktur, etwa aus einer Finite-Elemente-Diskretisierung (vgl. Seite 31).
171d.h. die Umnumerierung der Ecken, bei der die Bandbreite der Matrix minimal wird
172Das ist die mathematische Präzisierung für „wahrscheinlich unheimlich schwierig, auf jeden Fall genauso schwierig, wie alle schwierigsten Probleme, die man bisher gefunden hat“. Wir wollen das hier nicht weiter verfolgen. Wer mehr Informationen hierüber wünscht, lese nach bei M.R.GAREY und D.S.JOHNSON:
Computers and Intractibility", A Guide to the Theory of NP-Completeness, Freeman and Comp., San
−1 0 1 2 3 4 5 6 7
−2
−1 0 1 2 3 4
Diskretisierungsstruktur
Abbildung 83: Diskretisierungsstruktur
Um die zugehörige Diskretisierungsmatrix aufstellen zu können, benötigt man eine Nume-rierung der Knoten der Diskretisierung, bzw. der Ecken des Graphen174.
Wir wollen mit der folgenden (zugegebenerweise besonders dummen und fernliegenden) Numerierung beginnen:
−1 0 1 2 3 4 5 6 7
−2
−1 0 1 2 3 4
3 2 4
5 6
1 7
9 10
11
12 13
14
15
17 16
8
Abbildung 84: Erste Numerierung Dann hat die Matrix die folgende Struktur:
174Eine solche Numerierung ist ja tatsächlich schon Voraussetzung für die Dpeicherbarkeit der Daten.
0 5 10 15 0
2 4 6 8 10 12 14 16 18
nz = 79
Matrixstruktur zur Numerierung 1
Abbildung 85: Matrix zur ersten Numerierung
Wenden wir nun den Cuthill-McKee Algorithus an, werden wir zur folgenden Numerierung geführt
−1 0 1 2 3 4 5 6 7
−2
−1 0 1 2 3 4
Cuthill−McKee−Numerierung
1
2 3
4
5 6
7
8 9
10 11
12
13
14 15
16
17
Abbildung 86: Zweite Numerierung mit dem Ergebnis einer schon etwas kleineren Bandbreite.
0 5 10 15
0 2 4 6 8 10 12 14 16 18
nz = 81
Matrixstruktur zur ersten Cuthill−McKee−Ordnung
Abbildung 87: Matrix zur zweiten Numerierung
Wenn Sie die Numerierung einmal genau nachvollziehen, werden Sie bemerken, dass der
rung der Nachbarn von Ecke 1 durch den Algorithmus nicht festgelegt, welcher der beiden anschließen mit den Nummern 2 und 3 versehenen Ecken die 2 und welcher die 3 bekommen soll. Beide haben den Grad 4. In diesem Fall braucht man natürlich eine „tie-break-Regel“, die einen Fortgang des Verfahrens ermöglich. Über solche Regeln ist in den 70er Jahren (ohne große Fortschritte zu erzielen) viel diskutiert worden. Man kann z.B. die vormalige Ordnung in der letzten Numerierung verwenden.
Daß die endgültige Bandbreite stark von der Wahl des Anfangsknotens abhängt, sieht man etwa an dem folgenden zwei CM-Numerierungen:
0 1 2 3 4 5 6 7 8
−3
−2
−1 0 1 2 3
Graph mit zwei versch. CMK−Numerierungen
6 4 2 1 3 5 7
2 3 4 5 6 7
1
Abbildung 88: Zwei CM-Numerierungen
0 2 4 6 8
0 1 2 3 4 5 6 7 8
nz = 18
Matrix zur blauen Numerierung
0 2 4 6 8
0 1 2 3 4 5 6 7 8
nz = 18
Matrix zur roten Numerierung
Abbildung 89: Matrizen zu den beiden Numerierung
Dieses Beispiel legt sicher die Idee nahe, dass es nicht günstig ist, einen Startknoten zu wählen „der mitten im Graphen liegt“, sondern dass man vielmehr besser in einem „Rand-punkt“ des Graphen startet.
Tatsächlich strebt man als Startecke die Verwendung einer sogenanntenperipheren Ecke an. Dies ist eine Ecke mit maximaler Exzentrizität, wobei die Exzentrizität e(v) einer Ecke v definiert ist über
Definition der Exzentrizität e(v)
e(v) := max{d(v, w)|w ist Ecke vonG(A)}, und
d(v, w) =minimaler Abstand (in Kanten) der Ecken v und w.
Der Grund für dieses Ansinnen wird klar, wenn man beachtet, dass der Cuthill-McKee-Algorithmus die Numerierung in mehreren „Stufen“ vornimmt, deren Größe mit der Band-breite der zugehörigen Matrixdarstellung korelliert ist.
Definition 4.7 (Stufen-Struktur eines Graphen zur Wurzel v )
SeiGein (ungerichteter) Graph undM eine Teilmenge seiner Ecken. MitAdj(M) bezeich-nen wir dann die Menge der Ecken von G für die es verbindende Kanten in G zu einer Ecke aus M gibt und die selbst nicht zuM gehören.
Ist dann v eine Ecke von G, so definieren wir füri= 0,1,2, . . . die i-te Stufe zur Wurzelv, Si(v), induktiv über
S0(v) := {v},
S1(v) := Adj(S0(v)),
Si+1(v) := Adj(Si(v))\Si−1(v) für i= 1,2, . . . .
Mit L(v) bezeichnen wir den Index der letzten Stufe von v in G, die nicht leer ist, und nennen sie die Länge der Stufen-Struktur zu v.
Beispiel: In Abbildung 86 waren schon die Stufen zu Startknoten 1skizziert worden. Die mit durchbrochenen Linien umgebenen Ecken bilden die erste Stufe, die blau eingerahmten gehören der zweiten Stufe an und die nichteingerahmten sind die Mitglieder der letzten dritten Stufe. Es ist mithin L(1) = 3.
Beobachtungen:
(i) Wird der CM-Algorithmus mit der Ecke v gestartet, so numeriert er sukzessiv die Knoten der Stufen S0(v), S1(v), S2(v), . . . .
(ii) Es ist die Länge der Stufen-Struktur bei v gerade die Exzentrizität der Ecke v, L(v) =e(v).
(iii) Bezeichnen wir mit|Si(v)|die Anzahl der Elemente deri-ten Stufe vonGzur Wurzel v, so gilt für die Bandbreite der zur zugehörigen Numerierung gehörigen Matrix
m := max{|i−j| |aij ̸= 0} ≤ max
i:=1,...,e(v)(|Si−1(v)|+|Si(v)| −1). (255) Wegen der Beobachtung (iii) wird man die (maximale) Stufengröße möglichst klein machen wollen. Da die Gesamtzahl der Ecken auf die Stufen verteilt wird, erwartet man, dass dies der Fall sein wird, wenn die Länge der Stufenstruktur, L(v), möglichst groß ist175. Wegen (ii) möchte man deshalb eine Starteckev möglichst hoher Exzentrizität.
Nach (i) kann man die Stufen-Struktur und die Exzentrizität einer Ecke mit dem CM-Algorithmus bestimmen. Der folgende GPS-CM-Algorithmus zieht den CM-CM-Algorithmus dar-über hinaus zur Bestimmung einer „quasi-peripheren“ Ecke heran.
Algorithmus 4.8 (GPS-Algorithmus [Gibbs-Poole-Stockmeyer 76]) (1) Wähle eine beliebige Ecker von G(A).
(2) Starte den CM-Algorithmus mit Ecke r.
(3) Bestimme aus dem Ergebnis e(r).
(4) Wähle eine Ecke minimalen Gradesv in der letzten Stufe Se(r)(r).
(5) Starte CM mit v und bestimme e(v). Wenn e(v) > e(r) setze r := v und fahre bei (4) fort. Andernfalls
(6) Wähle v als quasi-peripheren Knoten und verwende ihn zur CM-Numerierung von G.
Zum einfacheren Verständnis dieses Algorithmus’ erlauben wir uns eine Reichlich unwissenschaftliche Interpretation des GPS-Algorithmus:
Fasse den Graphen als ein in den Ecken verknüpftes System von Lunten auf. Interpretiere die Stufen S0, S1, S2, . . . zu einer Startecke v als sukzessive Positionen einer Feuerfront zu den diskreten Zeitpunkten i= 0,1,2, . . ., die sich nach Anzünden des Systems in der Ecke v ausbreitet.
Ziel ist es, den „Zünd-Knoten“ ausfindig zu machen, bei dem das System am längsten brennt.
Iteration: Zünde das System irgendwo an. Stoppe die Brennzeit. Wähle aus der Stufe vor dem endgültigen Verlöschen die Ecke, die mit ihren Nachbarpunkten möglichst wenig Verbindungen hat, von der aus der Brand sich (also) am langsamsten auf das Gesamtsystem überträgt. Wiederhole damit das Vorgehen so lange, bis sich die Brennzeit nicht mehr verlängert.
Beispiel CM2: Sucht man einen Startpunkt für das System aus Beispiel CM2, so erhält man das folgende Ergebnis mit sechs Stufen (die wieder in das Bild eingezeichnet sind):
−1 0 1 2 3 4 5 6 7
−2
−1 0 1 2 3 4
Cuthill−McKee−Numerierung mit sechs Schichten
1 2
3 4
6 5 7
8 9 10
11 12
13
15 14
17 16
Abbildung 90: Dritte Numerierung
mit der anschließenden zugehörigen Matrixstruktur
0 5 10 15
0 2 4 6 8 10 12 14 16 18
nz = 81
Matrixstruktur zur zweiten Cuthill−McKee−Ordnung
Abbildung 91: Matrix zur dritten Numerierung Tatsächlich wird man die Numerierung am Ende noch umkehren
0 5 10 15
0 2 4 6 8 10 12 14 16 18
nz = 81
Zweite Cuthill−McKee−Ordnung, Revers
Abbildung 92: Matrix zur dritten Numerierung,rückwärts
weil man zeigen kann, dass sich dadurch die Größe der „Hülle der Matrix“ nicht vergrößert.
Die Hülle eine dünnbesetzten Matrix A = (aij)i,j=1,...,n erhält man als Verallgemeinerung ihres Bandes, indem man bei der Bandbreite „Zeilenabhängigkeit“ zulässt176.
Definition 4.9 (Bandbreite, Hülle)
Genauer definiert man zu A ∈ Rn×n, symmetrisch, und i ∈ {1, . . . , n} die (linke) Band-breite der i–ten Zeile durch
jM(i) := max{j |j ≤i, aik = 0∀k < j}. Dann ist die Hülle env(A) von A gegeben durch die Indexmenge
env(A) :={(i, j)|jM(i)≤j ≤i, i= 1, . . . , n}.
176In der anglo-amerikanischen Literatur heißen Hüllenmatrizen deshalb auch oft „variable bandwidth
Anmerkungen 4.10
Die Hülle ist die Menge an Matrix-Speicherplätzen, die durch die ursprünglichen Nicht-Null-Elemente belegt waren und für die man bei Gauß-Elimination ein Auffüllen nicht ausschließen kann.
Definition 4.11 (Fill-in)
Nicht-Null-Elemente, die während des Lösungsprozesses erzeugt werden, heißen „Fill-In“.
Achtung:Den allgemeinen Gepflogenheiten entsprechend wurde hier nur der linke untere Teil der Matrix betrachtet, weil man bei symmetrischen Matrizen üblicherweise auch nur diesen Teil speichert. Der Anteil oberhaupt der Hauptdiagonale ergibt sich durch Sym-metrie. Numerische Algorithmen der „Eliminationsklasse“ lassen sich zudem so schreiben, dass sie nur die Elemente der so „einseitig“ definierten Hülle benötigen.
Die nächste Figur zeigt links eine typische CM-Ergebnismatrix und rechts ihre durch so-genannten Fill-In komplettierte Hülle.
0 2 4 6 8 10
0 1 2 3 4 5 6 7 8 9 10
nz = 33
Typische Cuthill−McKee Struktur
0 2 4 6 8 10
0 1 2 3 4 5 6 7 8 9 10
nz = 33
Erzeugter "Fill−In"
Abbildung 93: CM-Matrix mit Fill-In
Eine Umkehr der Reihenfolge (RCM, „Reverse Cuthill-McKee“) der Numerierung wird die Hülle nachweislich nicht größer aber häufig wie im vorliegenden konstruierten Fall -deutlich kleiner machen:
0 2 4 6 8 10 0
1 2 3 4 5 6 7 8 9 10
nz = 33
Reverse Cuthill−McKee−Struktur
0 2 4 6 8 10
0 1 2 3 4 5 6 7 8 9 10
nz = 33
Fill−In bei Reverse CM
Abbildung 94: RCM-Matrix mit Fill-In