Cholesky-Zerlegung

Die Cholesky-Zerlegung (auch Cholesky-Faktorisierung) (nach André-Louis Cholesky, 1875–1918) bezeichnet in der numerischen Mathematik eine Zerlegung einer symmetrischen positiv definiten Matrix in ein Produkt aus einer unteren Dreiecksmatrix und deren Transponierter. Sie wurde von Cholesky vor 1914 im Zuge der Triangulation Kretas durch den Service géographique de l’armée entwickelt. Das Konzept kann auch allgemeiner für hermitesche Matrizen definiert werden.

Einsatzbereiche

Bei der Anwendung der Methode der kleinsten Quadrate ist eine Möglichkeit, die auftauchenden Minimierungsprobleme über die Normalgleichungen zu lösen, die eine symmetrische positiv definite Systemmatrix haben. Dies ist mit Hilfe der Cholesky-Zerlegung möglich und dies war die Motivation von Cholesky, die Zerlegung zu entwickeln. Beim Gauß-Newton-Verfahren ist damit bei jedem Iterationsschritt ein Gleichungssystem zu lösen, das sich mit dem Cholesky-Verfahren bestimmen lässt.

Die Cholesky-Zerlegung kann auch zur Gewinnung eines Vorkonditionierungsverfahrens für lineare Gleichungssysteme mit positiv definiter Matrix benutzt werden; zu diesem Zweck gibt es speziell die Varianten der unvollständigen Cholesky-Zerlegung und der modifizierten unvollständigen Cholesky-Zerlegung.

Gleichzeitig stellt die Zerlegung einen Test dar, ob eine gegebene symmetrische Matrix positiv definit ist. Andernfalls ist einer der Einträge auf der Hauptdiagonalen negativ, sodass die Wurzel nicht gezogen werden kann, oder gleich {\displaystyle 0}, sodass durch den Eintrag nicht dividiert werden kann. In beiden Fällen bricht der Algorithmus ab. Die Cholesky-Zerlegung lässt sich auch zur Bestimmung der Determinante der Matrix A verwenden, denn es gilt \det A=\prod _{{i=1}}^{n}G_{{ii}}^{2}.

Außerhalb der Mathematik findet die Cholesky-Zerlegung auch Anwendung in der ökonometrischen Erforschung makroökonomischer Zusammenhänge. Hierbei wird bei sogenannten vektorautoregressiven Modellen (VAR) die Reihenfolge der Beeinflussung der endogenen Variablen untereinander festgelegt.

Darüber hinaus wird sie auch bei der Monte-Carlo-Simulation eingesetzt, um vorgegebene Korrelationen in unabhängig generierte Zufallszahlenfolgen (als Diskretisierung stochastischer Prozesse) zu bringen.

Formulierung und Anwendung

Jede symmetrische, positiv definite Matrix A\in {\mathbb  {R}}^{{n\times n}} kann eindeutig in der Form

A=LDL^{{T}}

geschrieben werden. Dabei ist L> eine untere, normierte Dreiecksmatrix und D eine Diagonalmatrix mit positiven Einträgen. Mit der Quadratwurzel von D und dem Matrix-Faktor G, definiert durch

D=D^{{1/2}}D^{{1/2}}

und

G:=LD^{{1/2}},

wird die Cholesky-Zerlegung – äquivalent – auch formuliert als

A=GG^{{T}}.

Liegt eine Berechnung der Cholesky-Zerlegung vor, so lässt sich das Gleichungssystem Ax=b effizient durch Vorwärts- und Rückwärtseinsetzen lösen:

Berechnung

Setzt man {\displaystyle A=GG^{T}\in \mathbb {R} ^{n\times n}}, so erhält man für die Elemente von A=\left(a_{ij}\right)_{ij}:

{\displaystyle a_{ij}=\sum \limits _{k=1}^{j}g_{ik}g_{jk}\,,\quad i\geq j.}

Dieser Zusammenhang führt direkt auf die folgenden Formeln für {\displaystyle G=\left(g_{ij}\right)_{ij}}:

{\displaystyle g_{ij}={\begin{cases}0&\mathrm {f{\ddot {u}}r} \ i<j\\{\sqrt {a_{ii}-\sum \limits _{k=1}^{i-1}g_{ik}^{2}}}&\mathrm {f{\ddot {u}}r} \ i=j\\{\frac {1}{g_{jj}}}\left(a_{ij}-\sum \limits _{k=1}^{j-1}g_{ik}g_{jk}\right)&\mathrm {f{\ddot {u}}r} \ i>j.\end{cases}}}

Bei diesem Algorithmus ist es wichtig die Reihenfolge der berechneten Matrixeinträge richtig durchzuführen. Die Einträge werden spaltenweise berechnet und angefangen mit dem niedrigsten Zeilenindex.

Die Berechnung der Zerlegung {\displaystyle A=LDL^{T}} erfolgt in analoger Art und Weise für {\displaystyle L=\left(l_{ij}\right)_{ij}} und {\displaystyle D=\left(d_{ij}\right)_{ij}}:

{\displaystyle {\begin{aligned}d_{ij}&={\begin{cases}0&\mathrm {f{\ddot {u}}r} \ i\neq j\\a_{ii}-\sum _{k=1}^{i-1}l_{ik}^{2}d_{kk}&\mathrm {f{\ddot {u}}r} \ i=j\end{cases}}\\l_{ij}&={\begin{cases}0&\mathrm {f{\ddot {u}}r} \ i<j\\1&\mathrm {f{\ddot {u}}r} \ i=j\\{\frac {1}{d_{jj}}}\left(a_{ij}-\sum _{k=1}^{j-1}l_{ik}l_{jk}d_{kk}\right)&\mathrm {f{\ddot {u}}r} \;i>j.\end{cases}}\end{aligned}}}

Auch bei diesen Algorithmen ist es wichtig die Reihenfolge der berechneten Matrixeinträge richtig zu wählen. Zuerst muss man zum Index {\displaystyle j=1,\dotsc ,n} den Eintrag {\displaystyle d_{jj}} berechnen und anschließend die j-te Spalte der Matrix L, also: l_{{ij}} für {\displaystyle i=j+1,\dotsc ,n}.

Aufwand und Stabilität

Die Cholesky-Zerlegung ist numerisch stabil. Im Vergleich erfordert das Eliminationsverfahren nach Gauß mit seiner algorithmischen Umsetzung, der LR-Zerlegung, etwa doppelt so viele Operationen, da nicht nur eine Matrix G, sondern zwei Faktoren L und R berechnet werden müssen. Bei der Cholesky-Zerlegung treten {\displaystyle {\frac {1}{3}}n^{3}+O(n^{2})} arithmetische Operationen auf, davon {\tfrac  {1}{6}}n^{3} Multiplikationen, {\tfrac  {1}{2}}n^{2} Divisionen und n Wurzeloperationen.

Pseudocode

Die Berechnungen in obigen Formeln können in verschiedener Weise durchgeführt werden. Die nach Tadeusz Banachiewicz benannte Variante berechnet die untere Dreiecksmatrix zeilenweise. In Pseudocode sieht das Verfahren zur Zerlegung der Matrix A in die Form {\displaystyle GG^{T}} so aus:

Zugriffe (weiß) und Schreibvorgänge (gelb).
     For i = 1 To n
        For j = 1 To i
            Summe = a(i, j)
            For k = 1 To j-1
                Summe = Summe - a(i, k) * a(j, k)
            If i > j Then
                a(i, j) = Summe / a(j, j)   // Untere Dreiecksmatrix
            Else If Summe > 0 Then          // Diagonalelement
                a(i, i) = Sqrt(Summe)       // ... ist immer groesser Null
            Else
                ERROR                       // Die Matrix ist (wenigstens numerisch) nicht symmetrisch positiv definit

Die Laufindizes {\displaystyle i,j=1,\,,\ldots \,,n} im Code entsprechen der mathematischen Notierung von Elementen der Matrix {\displaystyle A=\left(a_{ij}\right)_{ij}}. Dabei ist n die Anzahl der Zeilen und gleichzeitig die Anzahl der Spalten der Matrix A, Hilfsvariablen sind k und Summe. Der Algorithmus arbeitet in-place; das heißt, er modifiziert die Matrix A so, dass sie die untere Dreiecksmatrix  G enthält. Es muss somit kein neuer Speicherplatz alloziert werden.

Der Algorithmus bearbeitet nur die linke untere Dreiecksmatrix von {\displaystyle A=\left(a_{ij}\right)_{ij}}, die Elemente  a_{ij} für i<j brauchen nicht mit Werten belegt zu werden (da die Matrix A nach Voraussetzung symmetrisch ist), und wenn sie Werte enthalten, werden diese nicht verändert. Sucht man also nach der Cholesky-Zerlegung  G gemäß {\displaystyle A=GG^{T}}, so sind die Matrixelemente von A oberhalb der Diagonalen (i<j) gleich 0 zu setzen.

Trenner
Basierend auf einem Artikel in: Wikipedia.de
Seitenende
Seite zurück
©  biancahoegel.de
Datum der letzten Änderung:  Jena, den: 04.04. 2023