Cholesky-Zerlegung

Die Cholesky-Zerlegung (nach André-Louis Cholesky (1875-1918)) bezeichnet in der numerischen Mathematik eine Zerlegung einer symmetrischen positiv definiten Matrix.

Formulierung und Anwendung

Jede symmetrische positiv definite Matrix $A\in \mathbb {R} ^{n\times n}$ ist orthogonal diagonalisierbar und kann daher eindeutig in der Form

$A=LDL^{T}\,$

geschrieben werden, wobei $L$ eine untere Dreiecksmatrix ohne Diagonale und $D$ eine positive Diagonalmatrix ist. Mit der Matrix-"Wurzel" $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 beispielsweise die Multiplikation der Matrixinversen $A^{-1}$ mit einem Vektor $b$ effizient durch Vorwärts- und Rückwärtseinsetzen berechnen:

• Durch Vorwärtseinsetzen Lösung des linearen Gleichungssystems $Gy=b$

• Durch Rückwärtseinsetzen Lösung des linearen Gleichungssystems $G^{T}x=y$

Damit kann die Choleskyzerlegung auch zur Gewinnung eines Vorkonditionierungsverfahrens für lineare Gleichungssysteme mit positiv definiter Matrix benutzt werden; zu diesem Zweck gibt es speziell die Variante der unvollständigen Cholesky-Zerlegung sowie der modifizierten unvollständigen Cholesky-Zerlegung. Es wird dabei ausgenutzt, dass man durch Cholesky-Zerlegung recht effizient eine (Näherungs-)Inverse der Matrix berechnen kann.

Berechnung

Formeln

Für alle k = 1..n; $g_{kk}={\sqrt {a_{kk}-\sum _{i=1}^{k-1}g_{ki}^{2}}}$

für l=k+1..n; Fehler beim Parsen (SVG (MathML kann über ein Browser-Plugin aktiviert werden): Ungültige Antwort („Math extension cannot connect to Restbase.“) von Server „http://localhost:6011/de.wikipedia.org/v1/“:): {\displaystyle g_{lk} = \frac{1}{g_{kk}} \cdot \left( a_{lk}-\sum_{i=1}^{k-1}g_{li} \cdot g_{ki} \right) }

Den Aufwand der Berechnung betreffend muß die Cholesky-Zerlegung mit dem Eliminationsverfahren nach Gauß und seiner algorithmischen Umsetzung, der LR-Zerlegung, verglichen werden. Letzteres erfordert etwa doppelt so viele Operationen, da eine Pivotisierung, d.h. Spalten- und/oder Zeilenpermutation, durchgeführt werden muß, einerseits um Division durch Null zu vermeiden, andererseits aus Stabilitätsgründen. Beide Gründe zur Pivotisierung entfallen bei positiv definiten Matrizen - dies ist die Ursache des geringeren Rechenaufwandes der Cholesky-Zerlegung. Die Cholesky-Zerlegung kann somit als eine Variante des Gauß-Verfahrens ohne Pivotisierung aufgefaßt werden.

Ein klassischer Testfall für jede Implementation der Cholesky-Zerlegung ist die Hilbert-Matrix, weil sie symmetrisch und positiv definit ist, aber eine vergleichsweise schlechte Kondition hat.

Bei klassischen Lösern linearer Gleichungssysteme stellt die Pivotisierung wiederum ein Hindernis für Algorithmen dar, welche die Bandbreite oder Hüllgröße des Faktors $G$ aus Gründen der Speichereffizienz minimieren sollen (etwa Cuthill-McKee- oder Minimum-Grad-Algorithmus). Demnach ist die Cholesky-Zerlegung auch besonders geeignet für die Zerlegung dünnbesetzter Matrizen.

Das Gelingen oder allfällige Nichtgelingen der Cholesky-Zerlegung ist ein wichtiges Kriterium um zu entscheiden, ob eine bestimmte Matrix positiv definit ist oder nicht.

Algorithmus

Die Berechnung eines linearen Gleichungssystems mit dem Cholesky-Verfahren kann in zwei Schritten vollzogen werden: Zunächst kann man die Matrix A (Dimension n × n), welche die Koeffizienten enthält, in die Form LDL^T bringen. In diesem Verfahrensschritt kann man feststellen, ob die Matrix A positiv und symmetrisch definit ist. Der Lösungsvektor B muss noch nicht bekannt sein.

In Pseudocode sieht das Cholesky-Verfahren zur Zerlegung der Matrix A in die Form LDL^T so aus:

   For i = 1 To n
       For k = i To n
           A(i,i) = A(i,i)-A(k,i)*A(k,i)
       Next k
       If A(i,i) <= 0 Then EXIT ("Matrix nicht Positiv definit!")
           else A(i,i)=sqrt(A(i,i))
       For j = i+1 To n
           For k = 1 To k-1
               A(i,j) = A(i,j) - A(k,i)*A(k,j)
           Next k
           A(i,j) = A(i,j) / A(i,i)
       Next j
   Next i

Die Indizes der Matrix A entsprechen der mathematischen Notierung i=1..n, j=1..n), n ist die Anzahl der Zeilen und gleichzeitig die Anzahl der Spalten der Matrix A, Hilfsvariablen sind i, j, k und Summe. Der Algorithmus arbeitet am Ort, d.h. er modifiziert die Matrix A so, dass sie die obere Dreiecksmatrix und die Diagonalmatrix DL^T enthält.

Der Algorithmus bearbeitet nur die rechte obere Dreiecksmatrix von A, die Werte $A_{i,j}\,$ für i < j brauchen nicht mit Werten belegt 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 $R$ gemäß $A=R^{T}R$, so sind die Matrixelemente von A unterhalb der Diagonalen (i < j) gleich Null zu setzen.

Sobald dann der Spaltenvektor B (Dimension n × 1) des Gleichungssystems zur Verfügung steht, kann weitergerechnet werden:

   For i = 1 To n
       For k = i - 1 To 1 Step -1
           b(i) = b(i) - a(i, k) * b(k)
       Next k
       b(i) = b(i) / a(i, i)
   Next i
   
   For i = n To 1 Step -1
       For k = i + 1 To n
           b(i) = b(i) - a(k, i) * b(k)
       Next k
       b(i) = b(i) / a(i, i)
   Next i

Der Vektor B wird durch die Berechnung mit dem Lösungsvektor überschrieben. Mit der vorhandenen Zerlegung von A können immer wieder Lösungen zu beliebig vielen verschiedenen Vektoren B bestimmt werden.

Einsatzbereiche

Bei der Anpassung von Polynomen oder Funktionen mit mehreren Veränderlichen, die mit der Gaußschen Fehlerquadratmethode vorgenommen werden, entstehen Gleichungssysteme, welche mit der Cholesky-Verfahren einfach gelöst werden können. Beim Gauß-Newton-Verfahren ist bei jedem Iterationsschritt ein Gleichungssystem zu lösen, welches sich mit dem Cholesky-Verfahren bestimmen lässt.

Vorgerechnetes Beispiel

$A=GG^{T}\,$

mit

$A={\begin{pmatrix}4&2&0\\2&5&2\\0&2&5\end{pmatrix}}$

und

$G={\begin{pmatrix}g_{11}&0&0\\g_{21}&g_{22}&0\\g_{31}&g_{32}&g_{33}\end{pmatrix}}$

und

$G^{T}={\begin{pmatrix}g_{11}&g_{21}&g_{31}\\0&g_{22}&g_{32}\\0&0&g_{33}\end{pmatrix}}$

ergibt sich

$G*G^{T}={\begin{pmatrix}g_{11}^{2}&g_{11}*g_{21}&g_{11}*g_{31}\\g_{21}*g_{11}&g_{21}^{2}+g_{22}^{2}&g_{21}*g_{31}+g_{22}*g_{32}\\g_{31}*g_{11}&g_{31}*g_{21}+g_{22}*g_{32}&g_{31}^{2}+g_{32}^{2}+g_{33}^{2}\end{pmatrix}}$.

Daraus folgt:

$g_{11}^{2}=a_{11}=4\rightarrow g_{11}=2$

$g_{21}*g_{11}=a_{21}=2\rightarrow g_{21}=1$

So wird für die restlichen g's weiterverfahren. Schließlich ergibt sich

$G={\begin{pmatrix}2&0&0\\1&2&0\\0&1&2\end{pmatrix}}$

und

$G^{T}={\begin{pmatrix}2&1&0\\0&2&1\\0&0&2\end{pmatrix}}$.

Literatur

• Hans Rudolf Schwarz & Nobert Köckler: Numerische Mathematik. Stuttgart: Teubner, 2004 (5. Auflage) ISBN 3-519-42960-8
• Gene H. Golub & Charles F. Van Loan: Matrix computations. Johns Hopkins University Press, 1996 (3rd edition) ISBN 0-80185414-8

Weblinks

Numerical Recipes in C Online Ausgabe