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