Factorizar la matriz $\mathbf{A}$
\begin{equation*} \mathbf{A} = \begin{bmatrix} 7 & 4 & 2 & 1 \\ 4 & 8 & 5 & 3 \\ 2 & 5 & 9 & 6 \\ 1 & 3 & 6 & 10 \end{bmatrix} \end{equation*}Calculando la primera columna
\begin{align*} l_{11} &= \sqrt{a_{11}} = \sqrt{7} = 2.645751 \\ l_{21} &= \frac{a_{21}}{l_{11}} = \frac{4}{2.645751} = 1.511858 \\ l_{31} &= \frac{a_{31}}{l_{11}} = \frac{2}{2.645751} = 0.755929 \\ l_{41} &= \frac{a_{41}}{l_{11}} = \frac{1}{2.645751} = 0.377964 \end{align*}Reemplazando
\begin{equation*} \begin{bmatrix} 2.645751 & 0 & 0 & 0 \\ 1.511858 & & 0 & 0 \\ 0.755929 & & & 0 \\ 0.377964 & & & \end{bmatrix} \end{equation*}Calculando la segunda columna
\begin{align*} l_{22} &= \sqrt{a_{22} - l_{21}^{2}} = \sqrt{8 - 1.511858^{2}} = 2.390457 \\ l_{32} &= \frac{a_{32} - l_{31} l_{21}}{l_{22}} = \frac{5 - 0.755929 (1.511858)}{2.390457} = 1.613559 \\ l_{42} &= \frac{a_{42} - l_{41} l_{21}}{l_{22}} = \frac{3 - 0.377964 (1.511858)}{2.390457} = 1.015945 \end{align*}Reemplazando
\begin{equation*} \begin{bmatrix} 2.645751 & 0 & 0 & 0 \\ 1.511858 & 2.390457 & 0 & 0 \\ 0.755929 & 1.613559 & & 0 \\ 0.377964 & 1.015945 & & \end{bmatrix} \end{equation*}Calculando la tercera columna
\begin{align*} l_{33} &= \sqrt{a_{33} - l_{31}^{2} - l_{32}^{2}} = \sqrt{9 - 0.755929^{2} - 1.613559^{2}} = 2.413503 \\ l_{43} &= \frac{a_{43} - l_{41} l_{31} - l_{42} l_{32}}{l_{33}} = \frac{6 - 0.377964 (0.755929) - 1.015945 (1.613559)}{2.413503} = 1.688417 \end{align*}Reemplazando
\begin{equation*} \begin{bmatrix} 2.645751 & 0 & 0 & 0 \\ 1.511858 & 2.390457 & 0 & 0 \\ 0.755929 & 1.613559 & 2.413503 & 0 \\ 0.377964 & 1.015945 & 1.688417 & \end{bmatrix} \end{equation*}Calculando la cuarta columna
\begin{equation*} l_{44} = \sqrt{a_{44} - l_{41}^{2} - l_{42}^{2} - l_{43}^{2}} = \sqrt{10 - 0.377964^{2} - 1.015945^{2} - 1.688417^{2}} = 2.444227 \end{equation*}Reemplazando
\begin{equation*} \begin{bmatrix} 2.645751 & 0 & 0 & 0 \\ 1.511858 & 2.390457 & 0 & 0 \\ 0.755929 & 1.613559 & 2.413503 & 0 \\ 0.377964 & 1.015945 & 1.688417 & 2.444227 \end{bmatrix} \end{equation*}Lo puede escribirse como
\begin{align*} l_{??} &= \sqrt{a_{??} - \sum_{k=1}^{?} l_{?k}^{2}} \\ l_{??} &= \frac{a_{??} - \sum_{k=1}^{?} l_{?k} l_{?k}}{l_{??}} \end{align*}para
\begin{align*} k &= 1 \\ &= 1, 2 \\ &= 1, 2, 3 \\ &= 1, 2, 3, 4 \end{align*}Lo puede escribirse como
\begin{align*} l_{??} &= \sqrt{a_{??} - \sum_{k=1}^{?} l_{?k}^{2}} \\ l_{i?} &= \frac{a_{i?} - \sum_{k=1}^{?} l_{ik} l_{?k}}{l_{??}} \end{align*}para
\begin{align*} i &= 2, 3, 4 \\ &= 3, 4 \\ &= 4 \end{align*}Lo puede escribirse como
\begin{align*} l_{jj} &= \sqrt{a_{jj} - \sum_{k=1}^{?} l_{jk}^{2}} \\ l_{ij} &= \frac{a_{ij} - \sum_{k=1}^{?} l_{ik} l_{jk}}{l_{jj}} \end{align*}para $j = 1, 2, 3, 4$
function cholesky_crout(a)
m, n = tamaño(a)
l = [[0,0,...,...,0],...,[0,0,...,...,0]]
for j=1 to n do
sumatoria = 0
for k=1 to j do
sumatoria = sumatoria + l(j,k)*l(j,k)
end for
l(j,j) = sqrt(a(j,j) - sumatoria)
for i=j+1 to m do
sumatoria = 0
for k=1 to j do
sumatoria = sumatoria + l(i,k)*l(j,k)
end for
l(i,j) = (a(i,j) - sumatoria)/l(j,j)
end for
end for
return l
end function
In [1]:
import numpy as np
def cholesky_crout(a):
m, n = a.shape
l = np.zeros((m,n))
for j in range(0,n):
sumatoria = 0
for k in range(0,j):
sumatoria = sumatoria + np.power(l[j,k],2)
l[j,j] = np.sqrt(a[j,j] - sumatoria)
for i in range(j+1,m):
sumatoria = 0
for k in range(0,j):
sumatoria = sumatoria + l[i,k]*l[j,k]
l[i,j] = (a[i,j] - sumatoria)/l[j,j]
return l
#versión alternativa
#def cholesky_crout(a):
# m, n = a.shape
# l = np.zeros((m,n))
# for j in range(0,n):
# l[j,j] = np.sqrt(a[j,j] - np.sum(np.power(l[j,:],2)))
# for i in range(j+1,m):
# l[i,j] = (a[i,j] - np.sum(l[i,:]*l[j,:]))/l[j,j]
# return l
In [2]:
A = np.array([[4,12,-16],[12,37,-43],[-16,-43,98]])
print(A)
In [3]:
cholesky_crout(A)
Out[3]:
In [4]:
#revisando el resultado
np.linalg.cholesky(A)
Out[4]:
In [5]:
B = np.array([[5,1.2,0.3,-0.6],
[1.2,6,-0.4,0.9],
[0.3,-0.4,8,1.7],
[-0.6,0.9,1.7,10]])
In [6]:
cholesky_crout(B)
Out[6]:
In [7]:
#revisando el resultado
np.linalg.cholesky(B)
Out[7]: