Algoritmo de Cholesky-Banachiewicz

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 fila

\begin{equation*} l_{11} = \sqrt{a_{11}} = \sqrt{7} = 2.645751 \end{equation*}

Reemplazando

\begin{equation*} \begin{bmatrix} 2.645751 & 0 & 0 & 0 \\ & & 0 & 0 \\ & & & 0 \\ & & & \end{bmatrix} \end{equation*}

Calculando la segunda fila

\begin{align*} l_{21} &= \frac{a_{21}}{l_{11}} = \frac{4}{2.645751} = 1.511858 \\ l_{22} &= \sqrt{a_{22} - l_{21}^{2}} = \sqrt{8 - 1.511858^{2}} = 2.390457 \end{align*}

Reemplazando

\begin{equation*} \begin{bmatrix} 2.645751 & 0 & 0 & 0 \\ 1.511858 & 2.390457 & 0 & 0 \\ & & & 0 \\ & & & \end{bmatrix} \end{equation*}

Calculando la tercera fila

\begin{align*} l_{31} &= \frac{a_{31}}{l_{11}} = \frac{2}{2.645751} = 0.755929 \\ l_{32} &= \frac{a_{32} - l_{31} l_{21}}{l_{22}} = \frac{5 - 0.755929 (1.511858)}{2.390457} = 1.613559 \\ l_{33} &= \sqrt{a_{33} - l_{31}^{2} - l_{32}^{2}} = \sqrt{9 - 0.755929^{2} - 1.613559^{2}} = 2.413503 \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 \\ & & & \end{bmatrix} \end{equation*}

Calculando la cuarta fila

\begin{align*} l_{41} &= \frac{a_{41}}{l_{11}} = \frac{1}{2.645751} = 0.377964 \\ l_{42} &= \frac{a_{42} - l_{41} l_{21}}{l_{22}} = \frac{3 - 0.377964 (1.511858)}{2.390457} = 1.015945 \\ 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 \\ 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{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 & 2.444227 \end{bmatrix} \end{equation*}

Patrón de cálculo

\begin{equation*} \begin{array}{llll} l_{11} = \sqrt{a_{11} - l_{11} l_{11}} & & & \\ l_{21} = \frac{a_{21} - l_{21} l_{11}}{l_{11}} & l_{22} = \sqrt{a_{22} - l_{21} l_{21} - l_{22} l_{22}} & & \\ 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} l_{31} - l_{32} l_{32} - l_{33} l_{33}} & \\ 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} l_{41} - l_{42} l_{42} - l_{43} l_{43} - l_{44} l_{44}} \end{array} \end{equation*}

Primer patrón

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

Lo puede escribirse como

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

Lo puede escribirse como

\begin{align*} l_{?j} &= \sqrt{a_{??} - \sum_{k=1}^{?} l_{?k} l_{jk}} \\ l_{?j} &= \frac{a_{?j} - \sum_{k=1}^{?} l_{?k} l_{jk}}{l_{jj}} \end{align*}

para

\begin{align*} j &= 1 \\ &= 1, 2 \\ &= 1, 2, 3 \\ &= 1, 2, 3, 4 \end{align*}

Tercer patrón

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

Lo puede escribirse como

\begin{align*} l_{ij} &= \sqrt{a_{ii} - \sum_{k=1}^{?} l_{ik} l_{jk}} \\ l_{ij} &= \frac{a_{ij} - \sum_{k=1}^{?} l_{ik} l_{jk}}{l_{jj}} \end{align*}

para $i = 1, 2, 3, 4 = 1, \dots, m$

Fórmula matemática

\begin{align*} i &= 1, \dots, m \\ & \quad j = 1, \dots, i \\ & \quad \quad l_{ij} = \begin{cases} \sqrt{a_{ii} - \sum_{k=1}^{j} l_{ik} l_{jk}} & \mbox{si } i = j \\ \cfrac{a_{ij} - \sum_{k=1}^{j} l_{ik} l_{jk}}{l_{jj}} & \mbox{si } i \neq j \end{cases} \end{align*}

Seudocódigo

function cholesky_banachiewicz(a)
    m, n = tamaño(a)
    l = [[0,0,...,...,0],...,[0,0,...,...,0]]
    for i=1 to m do
        for j=1 to i do
            sumatoria = 0
            for k=1 to j do
                sumatoria = sumatoria + l(i,k)*l(j,k)
            end for
        if(i=j) then
            l(i,j) = sqrt(a(i,i) - sumatoria)
        else
            l(i,j) = (a(i,j) - sumatoria)/l(j,j)
        end if
    end for
    return l
end function

Implementación



In [8]:

    
import numpy as np

def cholesky_banachiewicz(a):
    m, n = a.shape
    l = np.zeros((m,n))
    for i in range(0,m):
        for j in range(0,i+1):
            sumatoria = 0
            for k in range(0,j):
                sumatoria = sumatoria + l[i,k]*l[j,k]
            if i==j:
                l[i,j] = np.sqrt(a[i,i] - sumatoria)
            else:
                l[i,j] = (a[i,j] - sumatoria)/l[j,j]
    return l

#versión alternativa
#def cholesky_banachiewicz(a):
#    m, n = a.shape
#    l = np.zeros((m,n))
#    for i in range(0,m):
#        for j in range(0,i+1):
#            if i==j:
#                l[i,j] = np.sqrt(a[i,i] - np.sum(np.power(l[i,:j],2)))
#            else:
#                l[i,j] = (a[i,j] - np.sum(l[i,:j]*l[j,:j]))/l[j,j]
#    return l



In [9]:

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









    



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



In [10]:

    
cholesky_banachiewicz(A)









    Out[10]:





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



In [11]:

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









    Out[11]:





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



In [12]:

    
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 [13]:

    
cholesky_banachiewicz(B)









    Out[13]:





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



In [14]:

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









    Out[14]:





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