Algorithmus für die Cholesky-Zerlegung

Wir hatten die folgenden Zusammenhänge für die Einträge der Cholesky Zerlegung \(L\):

\[\begin{split}L_{kk} &= \sqrt{ A_{kk} - \sum_{j=1}^{k-1} L_{kj}^2 }\\ L_{ik} &= \frac 1{L_{kk}} \left(A_{kk} - \sum_{j=1}^{k-1} L_{kj}L_{ij} \right)\end{split}\]

Die führt uns auf die einfachste Variante eines Algorithmus:

def cholesky_classic(A):
    N = A.shape[0]
    L = numpy.zeros((N,N))
    for k in range(0,N):
         # compute diagonal entry
         L[k,k] = A[k,k]
         for j in range(0,k):
             L[k,k] = L[k,k] - L[k,j]*L[k,j]
         L[k,k] = sqrt(L[k,k])
         # compute remaining column
         for i in range(k+1,N):
             L[i,k] = A[i,k]
             for j in range(0,k):
                 L[i,k] = L[i,k] - L[i,j]*L[k,j]
             L[i,k] = L[i,k] / L[k,k]
    return L

Wir testen diese am vorherigen Beispiel:

print "classic"
A = numpy.array([[4., 2, 4, 4], [2, 10, 5, 2], [4, 5, 9, 6], [4, 2, 6, 9]])
print cholesky_classic(A)

Wenn wir uns die “händische” Konstruktion des vorangegangenen Beispiels vor Augen führen, so können wir die Einträge berechnen als:

def L_diag(A,k):
    L[k,k] = sqrt(A[k,k] - numpy.dot(L[k,0:k],L[k,0:k].transpose()))
def L_other(A,i,k):
    L[i,k] = 1.0/A[k,k] - (A[i,k] - numpy.dot(L[i,0:k],L[k,0:k].transpose()))

bitte beachten Sie wieder, dass Python die Indizierung von 0 bis N-1 führt und dass die Schreibweise 0:k dem Intervall \([0,k)\) entspricht.

Wir können dies zunächst nutzen um einen simplen (wenn auch nicht zwingend effizienten) Algorithmus zu konstruieren. Wir berechnen dafür zunächst die Skalarprodukte, dann ziehen wir die Wurzel des Diagonaleintrags und dann skalieren wir die ganze Spalte:

def cholesky_simple(A):
    N = A.shape[0]
    L = numpy.zeros((N,N))
    for k in range(0,N):
         L[k:,k] = A[k:,k] - numpy.dot(L[k:,0:k],L[k,0:k].transpose())
             # compute weight
         w = sqrt(L[k,k])
             # scale column
         L[k:,k] = 1.0/w * L[k:,k]
    return L

Wir testen die am letzten Beispiel:

print "simple"
print cholesky_simple(A)

Betrachten wir den Algorithmus genauer, so sehen wir, dass wir nur Einträge benötigen, die wir bereits berechnet haben, so dass die die Berechnung “inplace” durchführen können:

def cholesky_inplace(A):
    N = A.shape[0]
    for k in range(0,N):
         # compute sum
         A[k:,k] = A[k:,k] - numpy.dot(A[k:,0:k],A[k,0:k].transpose())
         # scale column
         w = sqrt(A[k,k])
         A[k:,k] = 1.0/w * A[k:,k]
         # clear upper triangle
         A[0:k,k] = 0
    return A

und erhalten das bekannte Ergebnis:

print "inplace"
A2 = numpy.copy(A)
print cholesky_inplace(A2)

Wir können den Algorithmus ausführlicher formulieren:

def cholesky_full(A):
    N = A.shape[0]
    for k in range(0,N):
         for i in range(k,N):
             for j in range(0,k):
                 A[i,k] = A[i,k] - A[i,j]*A[k,j]
         A[k,k] = sqrt(A[k,k])
         for i in range(k+1,N):
             A[i,k] = 1.0/A[k,k] * A[i,k]
         for i in range(0,k):
             A[i,k] = 0
    return A

print "full"
A2 = numpy.copy(A)
print cholesky_full(A2)

Im Speicher sind die Matrixeinträge jedoch meistens (z.B. in Python, C, C++, Java) zeilenweise gespeichert (in Fortran ist es spaltenweise). Daher ist es effizienter auch zeilenweise die Matrix \(L\) zu berechnen. Die folgende Variante berechnet (etwas weniger offensichtlich) die gleiche Zerlegung, allerdings findet die Berechnung zeilenweise statt:

def cholesky_rowwise(A):
    N = A.shape[0]
    for i in range(0,N):
        for k in range(0,i+1):
            sigma = A[i,k]
            for j in range(0,k):
                sigma = sigma - A[k,j] * A[i,j]
            if i > k:
                A[i,k] = sigma / A[k,k]
            else:
                A[k,k] = sqrt(sigma)
        for k in range(i+1,N):
                    A[i,k] = 0
    return A

print "row-wise"
A2 = numpy.copy(A)
print cholesky_rowwise(A2)

Dies ist der Algorithmus, wie wir ihn in der Vorlesung kennen gelernt hatten.

Und auch diesen Algorithmus kann man natürlich wieder mit Python Mitteln kompakter formulieren, allerdings nicht ganz so kompakt wie cholesky_inplace:

def cholesky_rowwise_compact(A):
    N = A.shape[0]
    for i in range(0,N):
        for k in range(0,i+1):
            A[i,k] = A[i,k] - numpy.dot(A[i,0:k],A[k,0:k])
            if i > k:
                A[i,k] = A[i,k] / A[k,k]
            else:
                A[k,k] = sqrt(A[k,k])
        A[i,i+1:N] = 0
    return A

print "row-wise compact"
A2 = numpy.copy(A)
print cholesky_rowwise_compact(A2)