Cholesky-Zerlegung Schritt für Schritt¶

Wir wollen in diesem Beispiel einmal selber nachvollziehen, was die Cholesky Zerlegung Schritt für Schritt berechnet. Wir wollen dazu die Zerlegung der folgenden \(4\times 4\) Matrix berechnen

\[\begin{split}A = \left( \begin{array}{cccc} 4 & 2 & 4 & 4 \\ 2 & 10 & 5 & 2 \\ 4 & 5 & 9 & 6 \\ 4 & 2 & 6 & 9 \end{array} \right)\end{split}\]

Wir legen die Matrix an:

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

und überprüfen, dass diese auch wirklich symmetrisch positiv definit ist, d.h. dass alle Eigenwerte größer 0 sind:

assert(numpy.all(numpy.linalg.eigvals(A) > 0))

dann legen wir eine leere Matrix \(L\) an:

L = numpy.zeros((4,4))

Wir erinnern uns, dass wir den Zusammenhang

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

Nach dieser Formel berechnen \(L_{11}\) (beachten Sie, dass Python die Indizierung bei 0 beginnt!):

print "compute diagonal entry 0"
L[0,0] = sqrt(A[0,0])
print L

und erhalten

\[\begin{split}\left( \begin{array}{cccc} \fbox{2} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{array} \right)\end{split}\]

Dann berechnen wir die restliche Spalte \(L_{21}\ldots L_{41}\):

# i = 1, k = 0
L[1,0] = 1.0/L[0,0] * A[1,0]
# i = 2, k = 0
L[2,0] = 1.0/L[0,0] * A[2,0]
# i = 3, k = 0
L[3,0] = 1.0/L[0,0] * A[3,0]

Wir können das auch einfacher schreiben als:

print "compute column 0"
L[1:4,0] = 1.0/L[0,0] * A[1:4,0]
print L

sodass direkt alle Einträge \(L_{21}\ldots L_{41}\) aufeinmal berechnet werden. Wir erhalten

\[\begin{split}\left( \begin{array}{cccc} 2 & 0 & 0 & 0 \\ \fbox{1} & 0 & 0 & 0 \\ \fbox{2} & 0 & 0 & 0 \\ \fbox{2} & 0 & 0 & 0 \end{array} \right)\end{split}\]

Nun können wir für \(k=2\) die zweite Spalte berechnen. Wir beginnen mit dem Diagonaleintrag:

print "compute diagonal entry 1"
L[1,1] = sqrt(A[1,1] - L[1,0]*L[1,0])
print L

\[\begin{split}\left( \begin{array}{cccc} 2 & 0 & 0 & 0 \\ 1 & \fbox{3} & 0 & 0 \\ 2 & 0 & 0 & 0 \\ 2 & 0 & 0 & 0 \end{array} \right)\end{split}\]

und berechnen anschließend wieder die restliche Spalte:

# i = 2, k = 1
L[2,1] = 1.0/L[1,1] * (A[2,1] - L[2,0]*L[1,0])
# i = 3, k = 1
L[3,1] = 1.0/L[1,1] * (A[3,1] - L[3,0]*L[1,0])

oder wieder in kurzer Schreibweise:

print "compute column 1"
L[2:4,1] = 1.0/L[1,1] * (A[2:4,1] - L[2:4,0]*L[1,0])
print L

\[\begin{split}\left( \begin{array}{cccc} 2 & 0 & 0 & 0 \\ 1 & 3 & 0 & 0 \\ 2 & \fbox{1} & 0 & 0 \\ 2 & \fbox{0} & 0 & 0 \end{array} \right)\end{split}\]

Analog berechnen wir zuerst \(L_{33}\) und dann den Rest der Spalte:

print "compute diagonal entry 2"
L[2,2] = sqrt(A[2,2] - (L[2,0]*L[2,0] + L[2,1]*L[2,1]))

wobei wir die Summe auch als Skalarprodukt auswerten können:

L[2,2] = sqrt(A[2,2] - numpy.dot(L[2,0:2],L[2,0:2].transpose()))
print L

\[\begin{split}\left( \begin{array}{cccc} 2 & 0 & 0 & 0 \\ 1 & 3 & 0 & 0 \\ 2 & 1 & \fbox{2} & 0 \\ 2 & 0 & 0 & 0 \end{array} \right)\end{split}\]

sowie \(L_{32}\):

L[3,2] = 1.0/L[2,2] * (A[3,2] - L[3,0]*L[2,0] - L[3,1]*L[2,1])

und wieder die Summe Skalarprodukt:

print "compute column 2"
L[3,2] = 1.0/L[2,2] * (A[3,2] - numpy.dot(L[3,0:1],L[2,0:1].transpose()))
print L

\[\begin{split}\left( \begin{array}{cccc} 2 & 0 & 0 & 0 \\ 1 & 3 & 0 & 0 \\ 2 & 1 & 2 & 0 \\ 2 & 0 & \fbox{1} & 0 \end{array} \right)\end{split}\]

Der letzte Eintrag ist wieder ein Diagonaleintrag:

print "compute diagonal entry 3"
L[3,3] = sqrt(A[3,3] - numpy.dot(L[3,0:3],L[3,0:3].transpose()))
print L

\[\begin{split}\left( \begin{array}{cccc} 2 & 0 & 0 & 0 \\ 1 & 3 & 0 & 0 \\ 2 & 1 & 2 & 0 \\ 2 & 0 & 1 & \fbox{2} \end{array} \right)\end{split}\]