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}$$