Wie berechnet man die Cholesky-Zerlegung?

Sei \(A \in \mathbb{R}^{n \times n}\) eine symmetrische, positiv definite Matrix. Dann existiert eine Zerlegung \(A = S \cdot D \cdot S^T\), wobei \(S\) eine unipotente Dreiecksmatrix ist und D eine positiv definite Diagonalmatrix.

Berechnung der Cholesky-Zerlegung

Hier ein paar Ausschnitte, aus der englischen Wikipedia:

Einfach von links oben nach rechts unten die Werte nach folgender Formel berechnen: \(D_j = A_{jj} - \sum_{k=1}^{j-1} S_{jk}^2 D_k\) \(S_{ij} = \frac{1}{D_j} \left( A_{ij} - \sum_{k=1}^{j-1} S_{ik} S_{jk} D_k \right), \qquad\text{for } i>j\)

Programmierer-Hinweise

Implementierung

Eine Python-Implementierung sieht so aus:

#!/usr/bin/env python
# -*- coding: utf-8 -*-

def getSD(A):
    """ @param A: eine quadratische, reele, positiv definite Matrix
        @return: Die Matrizen S und D, für die gilt:
                 A = S * D * S^T
    """
    n = len(A)
    S = [[0 for j in range(n)] for i in range(n)]
    D = [[0 for j in range(n)] for i in range(n)]

    for i in range(n):
        S[i][i] = 1.0

    for j in range(n):
        _summe  = sum(S[j][k]**2 * D[k][k] for k in range(j))
        D[j][j] = A[j][j] - _summe
        _summe  = sum(S[i][k]*S[j][k]*D[k][k] for k in range(j))
        S[i][j] = 1.0 / D[j][j] * (A[i][j] - _summe)
    return S, D

Bibliotheken

Ich habe mich mal nach Bibliotheken umgesehen, die die Cholesky-Zerlegung direkt beherrschen. NumPy kann es natürlich:

from numpy import linalg
print(linalg.cholesky([[5,1],[1,1]]))

Gibt aus:

array([[ 2.23606798,  0.        ],
       [ 0.4472136 ,  0.89442719]])

Das ist NICHT die Zerlegung \(A = S \cdot D \cdot S^T\), sondern \(A = G \cdot G^T\).

Interessanterweise hat NumPy das nicht direkt selbst implementiert (NumPy-Quelle). Man greift auf LAPACK zurück, was in FORTRAN 90 programmiert wurde (LAPACK-Quelle)!

Wolfram|Alpha

Auch Wolfram|Alpha kennt "cholesky decomposition": Beispiel. Allerdings sieht es schon bei \(3 \times 3\)-Matrizen schlecht aus.

Numerik

In Numerik haben wir bei Herrn Dr. Weiß folgendes als Cholesky-Zerlegung kennen gelernt:

Sei \(A \in \mathbb{R}^{n \times n}\) eine symmetrische, positiv definite Matrix. Dann existiert eine Zerlegung \(A = \bar L \cdot \bar{L}^T\), wobei \(\bar L\) eine untere Dreiecksmatrix ist.

Wenn man wie gewohnt eine LR-Zerlegung der Matrix \(A\) durchführt, erhält man zwei Matrizen \(L, R \in \mathbb{R}^{n \times n}\), wobei gilt: \(R = D \cdot L^T\), wobei \(D\) eine positiv definite Diagonalmatrix ist.

Offensichtlich gilt: \(\bar L = L \cdot D^{\frac{1}{2}}\).

Die Cholesky-Zerlegung kann man folgendermaßen berechnen:

In Python sieht das dann so aus:

def getL(A):
    n = len(A)
    L = [[0 for i in range(n)] for j in range(n)]
    print(L)
    print("")

    for k in range(n):
        L[k][k] = (A[k][k] - sum([L[k][i]**2 for i in range(k)]))**0.5
        for i in range(k+1, n):
            L[i][k] = (A[i][k] 
                        - sum([L[i][j]*L[k][j] for j in range(k)])) \
                      / L[k][k]
    return L

Siehe auch

• Cholesky decomposition (Englisch)
• Cholesky-Zerlegung