Das Lösen von linearen bzw. linearisierten Gleichungssystemen ist eine Aufgabe, die in fast allen Natur- und Ingenieurswissenschaften vorkommt. Aus diesem Grund existieren auch eine Reihe von nummerischen Bibliotheken wie bspw. das Basic Linear Algebra Subprograms (BLAS) oder das Linear Algebra PACKage (LAPACK), die für verschiedene Programmiersprachen auch portiert wurden. So existiert von LAPACK z.B. auch ein Pendant in JAVA unter dem Namen JLAPACK. Auch wenn der Einsatz dieser Bibliothek das Arbeiten mit Gleichungssystemen enorm vereinfacht, nehmen sie dem Anwender keine grundsätzlichen Überlegungen wie Speicherverbrauch oder Geschwindigkeit ab. In diesem Beitrag soll daher die Choleskey-Zerlegung zum Lösen von Gleichungssystemen näher betrachtet werden, da sie sich in vielfacher Hinsicht als vorteilhaft erweist (Knickmeyer 1987).

Die auf den gleichnamigen französischen Mathematiker zurückgehende Cholesky-Zerlegung zerlegt eine symmetrische positiv definite Matrix A in zwei Faktoren. Es gilt

Cholesky

mit

Zerlegung

Die Matrix L bzw. R entsprechen dabei der unteren (linken) bzw. oberen (rechten) Dreiecksmatrix. Das Lösen eines linearen bzw. linearisierten Gleichungssystems

Gleichung

wird in der (geodätischen) Fachliteratur häufig mittels der Inversen dargestellt.

Inverse

Das ist formell auch vollkommen korrekt, heißt aber noch lange nicht, dass es wirklich nötig ist, die Inverse von A explizit zu bestimmen um dieses System nach x aufzulösen. Das Programmsystem MATLAB (MATrix LABoratory) weist auf diesen Umstand in der Dokumentation auch explizit hin und verdeutlicht anhand eines Benchmarks den zeitlichen Mehraufwand beim Lösen mittels Inversion (Mathworks, www).

Die Cholesky-Zerlegung kann nur auf symmetrisch positiv definite Matrizen angewendet werden, was somit zunächst auch eine Einschränkung dieses Verfahrens ist. Berücksichtigt man jedoch, dass nahezu alle symmetrischen Matrizen im naturwissenschaftlich-technischen Bereich diese Eigenschaft erfüllen, kann dieser Umstand vernachlässigt werden (TM-Mathe, www). In der Geodäsie trifft man häufig im Zusammenhang mit homogenisierten Gleichungssystemen auf die Cholesky-Zerlegung. Hierbei wird das Gleichungssystem so transformiert, dass alle Beobachtungen unabhängig und gleichgenau sind (vgl. Caspary & Wichmann 2007). Für das Lösen von Gleichungssystemen bieten u.a. Ghilani & Wolf (2006) fertige Algorithmen in verschiedenen Programmiersprachen an. Darüber hinaus gehen die beiden Autoren auch auf die Speichertechniken für Matrizen ein und stellen einfache Konzepte vor. Des weiteren kann vollständiger C-Quellcode Engeln-Müllges & Reutter (1987) entnommen werden. Angelehnt an Knickmeyer (1987), der für große lineare Gleichungssysteme Berechnung der Inversen diskutiert, soll hier das Konzept zum Lösen von Gleichungssystemen und die algorithmische Ableitung erfolgen.

Das Lösen mittels Cholesky-Zerlegung gestaltet sich zweistufig und ist bei (TM-Mathe, www) anschaulich dargestellt. Der Einfachheit halber soll davon ausgegangen werden, dass die rechte, obere Dreiecksmatrix R aus der Zerlegung gewonnen wurde. Zunächst wird das Teilproblem

Vorwärtssub.

gelöst, welches anschließend durch

Rückwärtssub.

auf die gesuchte Lösung von x führt. Das Lösen dieser beiden Gleichungen wird häufig auch als Vorwärts- und Rückwärtssubstitution bezeichnet, da zum Auflösen nach y das System von oben nach unten gelöst wird und für die Bestimmung von x das Gleichungssystem von unten nach oben zu lösen ist. Ist auch die inverse Matrix von A gesucht, um bspw. Zuverlässigkeitsparameter abzuleiten, so kann diese abschließend aus

Inverse aus Zerlegung

gewonnen werden. Für die Bestimmung des Lösungsvektors ist sie jedoch nicht von Nöten.

Die Vorteile bei einer rechnergestützten Umsetzung erschließen sich ggf. noch nicht auf dem ersten Blick. Wird bedacht, dass die Matrix A eine symmetrische Matrix ist, genügt es, nur die untere oder obere Dreiecksmatrix zu speichern. Hierdurch gelingt somit eine Speicherreduzierung von fast 50%. Die Zerlegungsmatrix R ist wie oben bereits angedeutet ebenfalls nur eine Dreiecksmatrix. Da die Matrix A nach der Zerlegung nicht mehr benötigt wird, kann ein Algorithmus benutzt werden, der in-situ arbeitet. Dies bedeutet, die Matrix A wird bei der Cholesky-Zerlegung durch R überschrieben – es muss demnach kein weiterer Speicherplatz verfügbar gemacht werden. Da S, die inverse Matrix von R, ebenfalls eine Dreiecksmatrix ist, können bei deren Bestimmung ähnliche Überlegungen angestellt werden. Für sämtliche Berechnungen muss somit kein zusätzlicher Speicher verfügbar gemacht werden, sodass diese Art der Berechnung ressourcenschonend ist. Weitere Speicheroptimierungen können noch durch spezielle Umsortierungen erreicht werden (Caspary & Wichmann 2007).

Auch die Berechnungszeit profitiert natürlich von den halben Datensätzen, sodass die Cholesky-Zerlegung deutlich schneller als die LR-Zerlegung ist, bei der zwei Matrizen durch die Faktorisierung bestimmt werden. Durch die Dreiecksstruktur ergeben sich darüber hinaus einfache Rekursionsregeln, mit denen die Vorwärts- und Rückwärtssubstitution berechnet und letztlich die Parameter bestimmt werden. Diesen sollen im Folgenden kurz beschrieben werden.

Die Bestimmung der Elemente der Zerlegungsmatrix R erfolgt über

Dreiecksmatrix aus Cholesky-Zerlegung

Da sich die Lösung von y aus

Vorwärtssubstitutionslösung

ergibt und diese Gleichung identisch mit der dritten der Cholesky-Zerlegung ist, kann die Bestimmung von y mit der Faktorisierung zusammen ablaufen. Hierzu ist lediglich der Vektor b als neue letzte Spalte an die Matrix A zu schreiben. Die so erweiterte Matrix soll A‘ sein. Nach der Zerlegung ist die letzte Spalte in der erweiterten Matrix R‘ der Vektor y

Erweiterung der NGL

Damit die Rekursion beginnen kann, ist die Dreiecksstruktur der erweiterten Matrix wieder herzustellen. Die gewichtete Verbesserungsquadratsumme bei der Methode der kleinsten Quadrate lautet

Verbesserungsquadratsumme

und somit

Minimierungsfunktion

worin l der (reduzierte) Beobachtungsvektor und P die Gewichtsmatrix sind.

Da diese Gleichung identisch ist mit der zweiten Gleichung der Cholesky-Zerlegung, liefert die Zerlegung die Verbesserungsquadratsumme Ω wenn l^TPl als letztes Element in A‘ hinzugefügt wird.

Schematische Darstellung der erweiterten Matrix

Um abschließend die Lösung von x zu ermitteln, ist nachfolgende rekursive Auflösung nötig.

Bestimmung des Lösungsvektors

Soll neben x auch die Inverse von A berechnet werden, so ist zunächst S zu bestimmen. Wird S‘ dabei aus der erweiterten Matrix R‘ berechnet, so fällt der Vektor der Unbekannten x als letzte Spalte in S‘ wiederum direkt mit ab, wenn r_n+1,n+1 = 1 gesetzt wird.

Dreiecksmatrixerweiterung

Für die Invertierung von R bzw. der erweiterten Matrix R‘ lässt sich wiederum eine Rekursionsvorschrift finden, die leicht programmierbar ist.

Invertierungsvorschrift

Durch die abschließende Multiplikation von S mit sich selbst gelangt man zur Inversen von A.

Ableitung der Kovarianzmatrix aus Zerlegungsmatrix

Bei vielen nicht-linearen Aufgabenstellungen wird eine Lösung iterativ gewonnen, indem die ausgeglichenen Parameter x im nächsten Berechnungsschritt als neue Näherung eingeführt werden. Die Teillösungen der einzelnen durchgeführten Rechenschritte sind idR. uninteressant. Bei einer algorithmischen Umsetzung bietet es sich somit an, stets nur den Lösungsvektor x zu bestimmen und einzig im letzten Rechenschritt auch die inverse Matrix des Normalgleichungssystem für das Ableiten von Zuverlässigkeitsgrößen zu berechnen. Durch das Ausnutzen der speziellen Matrizenstrukturen arbeitet diese Variante nicht nur schnelle sondern auch ressourcenschonend.

Als Beispiel soll ein einfaches Höhennetz dienen, dass Caspary & Wichmann (2007) entnommen ist. Es handelt sich um ein zwangsfreies Netz mit einem Höhenanschluss (HB = 0.0000m).

Die Normalgleichung A und der Absolutgliedvektor b lauten für dieses Beispiel:

Normalgleichungssystem und Absolutgliedvektor

Die Zerlegung der erweiterten Matrix A‘ liefert zum einen die Verbesserungsquadratsumme Ω = 1.175 mm² und zum anderen die Teillösung y

Lösung der Vorwärtssubstitution

Um die Varianz-Kovarianz-Matrix für den Parametervektor x zu ermitteln, ist die Inverse von R zu ermitteln und mit sich selbst zu multiplizieren. Wie oben bereits dargestellt, kann die Berechnung von x und S über die erweiterte Matrix R‘ zusammen erfolgen. Die finale Lösung lautet

Lösung der Rückwärtssubstitution

5.0509731232623  -3.8461538461539                0                  0    -4.482224281742
              0   6.3461538461539             -2.5                  0   -16.053326923077
              0                 0  3.4708737864078   -0.9708737864078     2.928890776699
              0                 0                0    2.7890556045896    19.652086496028

5.05097     -3.84615      0.00000      0.00000      4.48222
0.00000      6.34615     -2.50000      0.00000     16.05333
0.00000      0.00000      3.47087     -0.97087     -2.92889
0.00000      0.00000      0.00000      2.78906    -19.65209
0.00000      0.00000      0.00000      0.00000    268.66062


2.24744     -1.71135      0.00000      0.00000      1.99437
0.00000      1.84863     -1.35235      0.00000     10.53018
0.00000      0.00000      1.28141     -0.75766      8.82748
0.00000      0.00000      0.00000      1.48829     -8.71058
0.00000      0.00000      0.00000      0.00000      0.00108


0.60560      0.53531      0.42717      0.14870     -7.13461
0.00000      0.70300      0.56098      0.19528     -8.20417
0.00000      0.00000      0.76684      0.26694     -3.42832
0.00000      0.00000      0.00000      0.45147      5.85274
0.00000      0.00000      0.00000      0.00000      1.00000

Die Cholesky-Zerlegung eignet sich aus mehreren Gründen zum Lösen von Gleichungssystemen. Die Zerlegung selbst lässt sich verhältnismäßig einfach programmieren. Durch das Ausnutzen der speziellen Matrizenstrukturen und in-situ Algorithmen ist kein zusätzlicher Speicher beim Lösen der Gleichungssysteme bereitzustellen. Gegenüber anderen Lösungsverfahren gilt die Cholesky-Zerlegung als nummerisch stabil und benötigt z.T. weniger Rechenoperationen als vergleichbare Lösungsansätze.

Verwendete Literatur, die nicht der Bibliothek entnommen ist:

• Dankert, J., Dankert, H. (www): Technische Mechanik - Cholesky-Verfahren, (zuletzt besucht: 14. Mai 2011).
• Knickmeyer, E.H.: Inversenberechung großer, linearer Gleichungssysteme. Deutsches Geodätisches Forschungsinstitut, München, 1979.
• MathWorks (www): Product Documentation: Matrix inverse (zuletzt besucht: 14. Mai 2011).

19.05.2011 von Michael Lösler