Algoritmo de Cholesky–Crout

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

Patrón de cálculo

$$ \begin{equation*} \begin{array}{llll} l_{11} = \sqrt{a_{11} - l_{11}^{2}} & & & \\ l_{21} = \frac{a_{21} - l_{21} l_{11}}{l_{11}} & l_{22} = \sqrt{a_{22} - l_{21}^{2} - l_{22}^{2}} & & \\ l_{31} = \frac{a_{31} - l_{31} l_{11}}{l_{11}} & l_{32} = \frac{a_{32} - l_{31} l_{21} - l_{32} l_{22}}{l_{22}} & l_{33} = \sqrt{a_{33} - l_{31}^{2} - l_{32}^{2} - l_{33}^{2}} & \\ l_{41} = \frac{a_{41} - l_{41} l_{11}}{l_{11}} & l_{42} = \frac{a_{42} - l_{41} l_{21} - l_{42} l_{22}}{l_{22}} & l_{43} = \frac{a_{43} - l_{41} l_{31} - l_{42} l_{32} - l_{43} l_{33}}{l_{33}} & l_{44} = \sqrt{a_{44} - l_{41}^{2} - l_{42}^{2} - l_{43}^{2} - l_{44}^{2}} \end{array} \end{equation*} $$

Primer patrón

$$ \begin{equation*} \begin{array}{llll} l_{11} = \sqrt{a_{11} - l_{1\color{blue}{1}}^{2}} & & & \\ l_{21} = \frac{a_{21} - l_{2\color{blue}{1}} l_{1\color{blue}{1}}}{l_{11}} & l_{22} = \sqrt{a_{22} - l_{2\color{blue}{1}}^{2} - l_{2\color{green}{2}}^{2}} & & \\ l_{31} = \frac{a_{31} - l_{3\color{blue}{1}} l_{1\color{blue}{1}}}{l_{11}} & l_{32} = \frac{a_{32} - l_{3\color{blue}{1}} l_{2\color{blue}{1}} - l_{3\color{green}{2}} l_{2\color{green}{2}}}{l_{22}} & l_{33} = \sqrt{a_{33} - l_{3\color{blue}{1}}^{2} - l_{3\color{green}{2}}^{2} - l_{3\color{red}{3}}^{2}} & \\ l_{41} = \frac{a_{41} - l_{4\color{blue}{1}} l_{1\color{blue}{1}}}{l_{11}} & l_{42} = \frac{a_{42} - l_{4\color{blue}{1}} l_{2\color{blue}{1}} - l_{4\color{green}{2}} l_{2\color{green}{2}}}{l_{22}} & l_{43} = \frac{a_{43} - l_{4\color{blue}{1}} l_{3\color{blue}{1}} - l_{4\color{green}{2}} l_{3\color{green}{2}} - l_{4\color{red}{3}} l_{3\color{red}{3}}}{l_{33}} & l_{44} = \sqrt{a_{44} - l_{4\color{blue}{1}}^{2} - l_{4\color{green}{2}}^{2} - l_{4\color{red}{3}}^{2} - l_{4\color{fuchsia}{4}}^{2}} \end{array} \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*} $$

Segundo patrón

$$ \begin{equation*} \begin{array}{llll} l_{11} = \sqrt{a_{11} - l_{11}^{2}} & & & \\ l_{\color{blue}{2}1} = \frac{a_{\color{blue}{2}1} - l_{\color{blue}{2}1} l_{11}}{l_{11}} & l_{22} = \sqrt{a_{22} - l_{21}^{2} - l_{22}^{2}} & & \\ l_{\color{green}{3}1} = \frac{a_{\color{green}{3}1} - l_{\color{green}{3}1} l_{11}}{l_{11}} & l_{\color{green}{3}2} = \frac{a_{\color{green}{3}2} - l_{\color{green}{3}1} l_{21} - l_{\color{green}{3}2} l_{22}}{l_{22}} & l_{33} = \sqrt{a_{33} - l_{31}^{2} - l_{32}^{2} - l_{33}^{2}} & \\ l_{\color{red}{4}1} = \frac{a_{\color{red}{4}1} - l_{\color{red}{4}1} l_{11}}{l_{11}} & l_{\color{red}{4}2} = \frac{a_{\color{red}{4}2} - l_{\color{red}{4}1} l_{21} - l_{\color{red}{4}2} l_{22}}{l_{22}} & l_{\color{red}{4}3} = \frac{a_{\color{red}{4}3} - l_{\color{red}{4}1} l_{31} - l_{\color{red}{4}2} l_{32} - l_{\color{red}{4}3} l_{33}}{l_{33}} & l_{44} = \sqrt{a_{44} - l_{41}^{2} - l_{42}^{2} - l_{43}^{2} - l_{44}^{2}} \end{array} \end{equation*} $$

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

Tercer patrón

$$ \begin{equation*} \begin{array}{llll} l_{\color{blue}{11}} = \sqrt{a_{\color{blue}{11}} - l_{\color{blue}{1}1}^{2}} & & & \\ l_{2\color{blue}{1}} = \frac{a_{2\color{blue}{1}} - l_{21} l_{\color{blue}{1}1}}{l_{\color{blue}{11}}} & l_{\color{green}{22}} = \sqrt{a_{\color{green}{22}} - l_{\color{green}{2}1}^{2} - l_{\color{green}{2}2}^{2}} & & \\ l_{3\color{blue}{1}} = \frac{a_{3\color{blue}{1}} - l_{31} l_{\color{blue}{1}1}}{l_{\color{blue}{11}}} & l_{3\color{green}{2}} = \frac{a_{3\color{green}{2}} - l_{31} l_{\color{green}{2}1} - l_{32} l_{\color{green}{2}2}}{l_{\color{green}{22}}} & l_{\color{red}{33}} = \sqrt{a_{\color{red}{33}} - l_{\color{red}{3}1}^{2} - l_{\color{red}{3}2}^{2} - l_{\color{red}{3}3}^{2}} & \\ l_{4\color{blue}{1}} = \frac{a_{4\color{blue}{1}} - l_{41} l_{\color{blue}{1}1}}{l_{\color{blue}{11}}} & l_{4\color{green}{2}} = \frac{a_{4\color{green}{2}} - l_{41} l_{\color{green}{2}1} - l_{42} l_{\color{green}{2}2}}{l_{\color{green}{22}}} & l_{4\color{red}{3}} = \frac{a_{4\color{red}{3}} - l_{41} l_{\color{red}{3}1} - l_{42} l_{\color{red}{3}2} - l_{43} l_{\color{red}{3}3}}{l_{\color{red}{33}}} & l_{\color{fuchsia}{44}} = \sqrt{a_{\color{fuchsia}{44}} - l_{\color{fuchsia}{4}1}^{2} - l_{\color{fuchsia}{4}2}^{2} - l_{\color{fuchsia}{4}3}^{2} - l_{\color{fuchsia}{4}4}^{2}} \end{array} \end{equation*} $$

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 $

Fórmula matemática

$$ \begin{align*} j &= 1, \dots, n \\ & \quad l_{jj} = \sqrt{a_{jj} - \sum_{k=1}^{j} l_{jk}^{2}} \\ & \quad i = 1 + j, \dots, m \\ & \quad \quad l_{ij} = \frac{a_{ij} - \sum_{k=1}^{j} l_{ik} l_{jk}}{l_{jj}} \end{align*} $$

Seudocódigo

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

Implementación

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

A = np.array([[4,12,-16],[12,37,-43],[-16,-43,98]])
print(A)

    [[  4  12 -16]
     [ 12  37 -43]
     [-16 -43  98]]

cholesky_crout(A)

    array([[ 2.,  0.,  0.],
           [ 6.,  1.,  0.],
           [-8.,  5.,  3.]])

#revisando el resultado
np.linalg.cholesky(A)

    array([[ 2.,  0.,  0.],
           [ 6.,  1.,  0.],
           [-8.,  5.,  3.]])

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]])

cholesky_crout(B)

    array([[ 2.23606798,  0.        ,  0.        ,  0.        ],
           [ 0.53665631,  2.38997908,  0.        ,  0.        ],
           [ 0.13416408, -0.19749127,  2.81833234,  0.        ],
           [-0.26832816,  0.43682391,  0.64657701,  3.05272387]])

#revisando el resultado
np.linalg.cholesky(B)

    array([[ 2.23606798,  0.        ,  0.        ,  0.        ],
           [ 0.53665631,  2.38997908,  0.        ,  0.        ],
           [ 0.13416408, -0.19749127,  2.81833234,  0.        ],
           [-0.26832816,  0.43682391,  0.64657701,  3.05272387]])