Algoritmo de Crout

$$ \begin{equation*} \begin{bmatrix} l_{11} & 0 & 0 & 0 \\ l_{21} & l_{22} & 0 & 0 \\ l_{31} & l_{32} & l_{33} & 0 \\ l_{41} & l_{42} & l_{43} & l_{44} \end{bmatrix} \begin{bmatrix} 1 & u_{12} & u_{13} & u_{14} \\ 0 & 1 & u_{23} & u_{24} \\ 0 & 0 & 1 & u_{34} \\ 0 & 0 & 0 & 1 \end{bmatrix} = \begin{bmatrix} a_{11} & a_{12} & a_{13} & a_{14} \\ a_{21} & a_{22} & a_{23} & a_{24} \\ a_{31} & a_{32} & a_{33} & a_{34} \\ a_{41} & a_{42} & a_{43} & a_{44} \end{bmatrix} \end{equation*} $$

Multiplicando $$ \begin{equation*} \begin{bmatrix} l_{11} & l_{11} u_{12} & l_{11} u_{13} & l_{11} u_{14} \\ l_{21} & l_{21} u_{12} + l_{22} & l_{21} u_{13} + l_{22} u_{23} & l_{21} u_{14} + l_{22} u_{24} \\ l_{31} & l_{31} u_{12} + l_{32} & l_{31} u_{13} + l_{32} u_{23} + l_{33} & l_{31} u_{14} + l_{32} u_{24} + l_{33} u_{34} \\ l_{41} & l_{41} u_{12} + l_{42} & l_{41} u_{13} + l_{42} u_{23} + l_{43} & l_{41} u_{14} + l_{42} u_{24} + l_{43} u_{34} + l_{44} \end{bmatrix} = \begin{bmatrix} a_{11} & a_{12} & a_{13} & a_{14} \\ a_{21} & a_{22} & a_{23} & a_{24} \\ a_{31} & a_{32} & a_{33} & a_{34} \\ a_{41} & a_{42} & a_{43} & a_{44} \end{bmatrix} \end{equation*} $$

Despejando y reemplazando $$ \begin{equation} \begin{bmatrix} a_{11} & 0 & 0 & 0 \\ a_{21} & a_{22} - l_{21} u_{12} & 0 & 0 \\ a_{31} & a_{32} - l_{31} u_{12} & a_{33} - l_{31} u_{13} - l_{32} u_{23} & 0 \\ a_{41} & a_{42} - l_{41} u_{12} & a_{43} - l_{41} u_{13} - l_{42} u_{23} & a_{44} - l_{41} u_{14} - l_{42} u_{24} - l_{43} u_{34} \end{bmatrix} \begin{bmatrix} 1 & \frac{a_{12}}{l_{11}} & \frac{a_{13}}{l_{11}} & \frac{a_{14}}{l_{11}} \\ 0 & 1 & \frac{a_{23} - l_{21} u_{13}}{l_{22}} & \frac{a_{24} - l_{21} u_{14}}{l_{22}} \\ 0 & 0 & 1 & \frac{a_{34} - l_{31} u_{14} - l_{32} u_{24}}{l_{33}} \\ 0 & 0 & 0 & 1 \end{bmatrix} = \begin{bmatrix} a_{11} & a_{12} & a_{13} & a_{14} \\ a_{21} & a_{22} & a_{23} & a_{24} \\ a_{31} & a_{32} & a_{33} & a_{34} \\ a_{41} & a_{42} & a_{43} & a_{44} \end{bmatrix} \tag{8.1} \end{equation} $$

Ejemplo

Factorizar la matriz $ \mathbf{A} $ $$ \begin{equation*} \begin{bmatrix} 1 & 1 & 2 & 3 \\ 2 & 1 & -1 & 1 \\ 3 & -1 & -1 & 2 \\ -1 & 2 & 3 & -1 \end{bmatrix} \end{equation*} $$

Calculando la primera columna de $ \mathbf{L} $ y la primera fila $ \mathbf{U} $ $$ \begin{align*} l_{11} &= a_{11} = 1 \\ l_{21} &= a_{21} = 2 \\ l_{31} &= a_{31} = 3 \\ l_{41} &= a_{41} = -1 \\ u_{12} &= \frac{a_{12}}{l_{11}} = \frac{1}{1} = 1 \\ u_{13} &= \frac{a_{13}}{l_{11}} = \frac{2}{1} = 2 \\ u_{14} &= \frac{a_{14}}{l_{11}} = \frac{3}{1} = 3 \end{align*} $$

Reemplazando $$ \begin{equation*} \begin{bmatrix} 1 & 0 & 0 & 0 \\ 2 & & 0 & 0 \\ 3 & & & 0 \\ -1 & & & \end{bmatrix} \begin{bmatrix} 1 & 1 & 2 & 3 \\ 0 & 1 & & \\ 0 & 0 & 1 & \\ 0 & 0 & 0 & 1 \end{bmatrix} \end{equation*} $$

Calculando la segunda columna de $ \mathbf{L} $ y la segunda fila $ \mathbf{U} $ $$ \begin{align*} l_{22} &= a_{22} - l_{21} u_{12} = 1 - 2 (1) = -1 \\ l_{32} &= a_{32} - l_{31} u_{12} = -1 - 3 (1) = -4 \\ l_{42} &= a_{42} - l_{41} u_{12} = 2 - (-1) (1) = 3 \\ u_{23} &= \frac{a_{23} - l_{21} u_{13}}{l_{22}} = \frac{-1 - 2 (2)}{-1} = 5 \\ u_{24} &= \frac{a_{24} - l_{21} u_{14}}{l_{22}} = \frac{1 - 2 (3)}{-1} = 5 \end{align*} $$

Reemplazando $$ \begin{equation*} \begin{bmatrix} 1 & 0 & 0 & 0 \\ 2 & -1 & 0 & 0 \\ 3 & -4 & & 0 \\ -1 & 3 & & \end{bmatrix} \begin{bmatrix} 1 & 1 & 2 & 3 \\ 0 & 1 & 5 & 5 \\ 0 & 0 & 1 & \\ 0 & 0 & 0 & 1 \end{bmatrix} \end{equation*} $$

Calculando la tercera columna de $ \mathbf{L} $ y la tercera fila $ \mathbf{U} $ $$ \begin{align*} l_{33} &= a_{33} - l_{31} u_{13} - l_{32} u_{23} = -1 - 3 (2) - (-4) (5) = 13 \\ l_{43} &= a_{43} - l_{41} u_{13} - l_{42} u_{23} = 3 - (-1) (2) - 3 (5) = -10 \\ u_{34} &= \frac{a_{34} - l_{31} u_{14} - l_{32} u_{24}}{l_{33}} = \frac{2 - 3 (3) - (-4) (5)}{13} = 1 \end{align*} $$

Reemplazando $$ \begin{equation*} \begin{bmatrix} 1 & 0 & 0 & 0 \\ 2 & -1 & 0 & 0 \\ 3 & -4 & 13 & 0 \\ -1 & 3 & -10 & \end{bmatrix} \begin{bmatrix} 1 & 1 & 2 & 3 \\ 0 & 1 & 5 & 5 \\ 0 & 0 & 1 & 1 \\ 0 & 0 & 0 & 1 \end{bmatrix} \end{equation*} $$

Calculando la cuarta columna de $ \mathbf{L} $ $$ \begin{equation*} l_{44} = a_{44} - l_{41} u_{14} - l_{42} u_{24} - l_{43} u_{34} = -1 - (-1) (3) - 3 (5) - (-10) (1) = -3 \end{equation*} $$

Reemplazando $$ \begin{equation*} \begin{bmatrix} 1 & 0 & 0 & 0 \\ 2 & -1 & 0 & 0 \\ 3 & -4 & 13 & 0 \\ -1 & 3 & -10 & -3 \end{bmatrix} \begin{bmatrix} 1 & 1 & 2 & 3 \\ 0 & 1 & 5 & 5 \\ 0 & 0 & 1 & 1 \\ 0 & 0 & 0 & 1 \end{bmatrix} \end{equation*} $$

Patrón de cálculo

Primer patrón

$$ \begin{equation*} \begin{array}{llll} \color{blue}{l_{11} = a_{11}} & u_{12} = \cfrac{a_{12}}{l_{11}} & u_{13} = \cfrac{a_{13}}{l_{11}} & u_{14} = \cfrac{a_{14}}{l_{11}} \\ \color{blue}{l_{21} = a_{21}} & l_{22} = a_{22} - l_{21} u_{12} & u_{23} = \cfrac{a_{23} - l_{21} u_{13}}{l_{22}} & u_{24} = \cfrac{a_{24} - l_{21} u_{14}}{l_{22}} \\ \color{blue}{l_{31} = a_{31}} & l_{32} = a_{32} - l_{31} u_{12} & l_{33} = a_{33} - l_{31} u_{13} - l_{32} u_{23} & u_{34} = \cfrac{a_{34} - l_{31} u_{14} - l_{32} u_{24}}{l_{33}} \\ \color{blue}{l_{41} = a_{41}} & l_{42} = a_{42} - l_{41} u_{12} & l_{43} = a_{43} - l_{41} u_{13} - l_{42} u_{23} & l_{44} = a_{44} - l_{41} u_{14} - l_{42} u_{24} - l_{43} u_{34} \end{array} \end{equation*} $$

Lo anterior puede escribirse como $$ \begin{equation*} l_{i1} = a_{i1} \end{equation*} $$

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

Segundo patrón

$$ \begin{equation*} \begin{array}{llll} l_{11} = a_{11} & \color{blue}{u_{12} = \cfrac{a_{12}}{l_{11}}} & \color{blue}{u_{13} = \cfrac{a_{13}}{l_{11}}} & \color{blue}{u_{14} = \cfrac{a_{14}}{l_{11}}} \\ l_{21} = a_{21} & l_{22} = a_{22} - l_{21} u_{12} & u_{23} = \cfrac{a_{23} - l_{21} u_{13}}{l_{22}} & u_{24} = \cfrac{a_{24} - l_{21} u_{14}}{l_{22}} \\ l_{31} = a_{31} & l_{32} = a_{32} - l_{31} u_{12} & l_{33} = a_{33} - l_{31} u_{13} - l_{32} u_{23} & u_{34} = \cfrac{a_{34} - l_{31} u_{14} - l_{32} u_{24}}{l_{33}} \\ l_{41} = a_{41} & l_{42} = a_{42} - l_{41} u_{12} & l_{43} = a_{43} - l_{41} u_{13} - l_{42} u_{23} & l_{44} = a_{44} - l_{41} u_{14} - l_{42} u_{24} - l_{43} u_{34} \end{array} \end{equation*} $$

Lo anterior puede escribirse como $$ \begin{equation*} l_{1j} = \frac{a_{1j}}{l_{11}} \end{equation*} $$

para $ j=2, 3, 4 = 2, \dots, n $

Tercer patrón

$$ \begin{equation*} \begin{array}{llll} l_{11} = a_{11} & u_{12} = \cfrac{a_{12}}{l_{11}} & u_{13} = \cfrac{a_{13}}{l_{11}} & u_{14} = \cfrac{a_{14}}{l_{11}} \\ l_{21} = a_{21} & l_{\color{blue}{22}} = a_{\color{blue}{22}} - l_{\color{blue}{2}1} u_{1\color{blue}{2}} & u_{\color{blue}{2}3} = \cfrac{a_{\color{blue}{2}3} - l_{\color{blue}{2}1} u_{13}}{l_{\color{blue}{22}}} & u_{\color{blue}{2}4} = \cfrac{a_{\color{blue}{2}4} - l_{\color{blue}{2}1} u_{14}}{l_{\color{blue}{22}}} \\ l_{31} = a_{31} & l_{3\color{blue}{2}} = a_{3\color{blue}{2}} - l_{31} u_{1\color{blue}{2}} & l_{\color{green}{33}} = a_{\color{green}{33}} - l_{\color{green}{3}1} u_{1\color{green}{3}} - l_{\color{green}{3}2} u_{2\color{green}{3}} & u_{\color{green}{3}4} = \cfrac{a_{\color{green}{3}4} - l_{\color{green}{3}1} u_{14} - l_{\color{green}{3}2} u_{24}}{l_{\color{green}{33}}} \\ l_{41} = a_{41} & l_{4\color{blue}{2}} = a_{4\color{blue}{2}} - l_{41} u_{1\color{blue}{2}} & l_{4\color{green}{3}} = a_{4\color{green}{3}} - l_{41} u_{1\color{green}{3}} - l_{42} u_{2\color{green}{3}} & l_{44} = a_{44} - l_{41} u_{14} - l_{42} u_{24} - l_{43} u_{34} \end{array} \end{equation*} $$

Lo anterior puede escribirse como $$ \begin{align*} l_{jj} &= a_{jj} - \sum l_{j?} u_{?j} \\ l_{?j} &= a_{?j} - \sum l_{??} u_{?j} \\ u_{j?} &= \frac{a_{j?} - \sum l_{j?} u_{??}}{l_{jj}} \end{align*} $$

para $ j=2, 3 = 2, \dots, n-1 $

Cuarto patrón

$$ \begin{equation*} \begin{array}{llll} l_{11} = a_{11} & u_{12} = \cfrac{a_{12}}{l_{11}} & u_{13} = \cfrac{a_{13}}{l_{11}} & u_{14} = \cfrac{a_{14}}{l_{11}} \\ l_{21} = a_{21} & l_{22} = a_{22} - l_{21} u_{12} & u_{2\color{blue}{3}} = \cfrac{a_{2\color{blue}{3}} - l_{21} u_{1\color{blue}{3}}}{l_{22}} & u_{2\color{blue}{4}} = \cfrac{a_{2\color{blue}{4}} - l_{21} u_{1\color{blue}{4}}}{l_{22}} \\ l_{31} = a_{31} & l_{\color{blue}{3}2} = a_{\color{blue}{3}2} - l_{\color{blue}{3}1} u_{12} & l_{33} = a_{33} - l_{31} u_{13} - l_{32} u_{23} & u_{3\color{green}{4}} = \cfrac{a_{3\color{green}{4}} - l_{31} u_{1\color{green}{4}} - l_{32} u_{2\color{green}{4}}}{l_{33}} \\ l_{41} = a_{41} & l_{\color{blue}{4}2} = a_{\color{blue}{4}2} - l_{\color{blue}{4}1} u_{12} & l_{\color{green}{4}3} = a_{\color{green}{4}3} - l_{\color{green}{4}1} u_{13} - l_{\color{green}{4}2} u_{23} & l_{44} = a_{44} - l_{41} u_{14} - l_{42} u_{24} - l_{43} u_{34} \end{array} \end{equation*} $$

Lo anterior puede escribirse como $$ \begin{align*} l_{ij} &= a_{ij} - \sum l_{i?} u_{?j} \\ u_{ji} &= \frac{a_{ji} - \sum l_{j?} u_{?i}}{l_{jj}} \end{align*} $$

para $$ \begin{align*} i &= 3, 4 \\ &= 4 \end{align*} $$

Quinto patrón

$$ \begin{equation*} \begin{array}{llll} l_{11} = a_{11} & u_{12} = \cfrac{a_{12}}{l_{11}} & u_{13} = \cfrac{a_{13}}{l_{11}} & u_{14} = \cfrac{a_{14}}{l_{11}} \\ l_{21} = a_{21} & l_{22} = a_{22} - l_{2\color{blue}{1}} u_{\color{blue}{1}2} & u_{23} = \cfrac{a_{23} - l_{2\color{blue}{1}} u_{\color{blue}{1}3}}{l_{22}} & u_{24} = \cfrac{a_{24} - l_{2\color{blue}{1}} u_{\color{blue}{1}4}}{l_{22}} \\ l_{31} = a_{31} & l_{32} = a_{32} - l_{3\color{blue}{1}} u_{\color{blue}{1}2} & l_{33} = a_{33} - l_{3\color{green}{1}} u_{\color{green}{1}3} - l_{3\color{green}{2}} u_{\color{green}{2}3} & u_{34} = \cfrac{a_{34} - l_{3\color{green}{1}} u_{\color{green}{1}4} - l_{3\color{green}{2}} u_{\color{green}{2}4}}{l_{33}} \\ l_{41} = a_{41} & l_{42} = a_{42} - l_{4\color{blue}{1}} u_{\color{blue}{1}2} & l_{43} = a_{43} - l_{4\color{green}{1}} u_{\color{green}{1}3} - l_{4\color{green}{2}} u_{\color{green}{2}3} & l_{44} = a_{44} - l_{41} u_{14} - l_{42} u_{24} - l_{43} u_{34} \end{array} \end{equation*} $$

Lo anterior puede escribirse como $$ \begin{align*} l_{jj} &= a_{jj} - \sum l_{jk} u_{kj} \\ l_{ij} &= a_{ij} - \sum l_{ik} u_{kj} \\ u_{ji} &= \frac{a_{ji} - \sum l_{jk} u_{ki}}{l_{jj}} \end{align*} $$

para $$ \begin{align*} k &= 1 \\ &= 1, 2 \end{align*} $$

Sexto patrón

$$ \begin{equation*} \begin{array}{llll} l_{11} = a_{11} & u_{12} = \cfrac{a_{12}}{l_{11}} & u_{13} = \cfrac{a_{13}}{l_{11}} & u_{14} = \cfrac{a_{14}}{l_{11}} \\ l_{21} = a_{21} & l_{22} = a_{22} - l_{21} u_{12} & u_{23} = \cfrac{a_{23} - l_{21} u_{13}}{l_{22}} & u_{24} = \cfrac{a_{24} - l_{21} u_{14}}{l_{22}} \\ l_{31} = a_{31} & l_{32} = a_{32} - l_{31} u_{12} & l_{33} = a_{33} - l_{31} u_{13} - l_{32} u_{23} & u_{34} = \cfrac{a_{34} - l_{31} u_{14} - l_{32} u_{24}}{l_{33}} \\ l_{41} = a_{41} & l_{42} = a_{42} - l_{41} u_{12} & l_{43} = a_{43} - l_{41} u_{13} - l_{42} u_{23} & \color{blue}{l_{44} = a_{44} - l_{41} u_{14} - l_{42} u_{24} - l_{43} u_{34}} \end{array} \end{equation*} $$

Lo anterior puede escribirse como $$ \begin{equation*} l_{nn} = a_{nn} - \sum_{k=1}^{n-1} l_{nk} u_{kn} \end{equation*} $$

Fórmula matemática

$$ \begin{align*} i &= 1, \dots, m \\ & \quad l_{i1} = a_{i1} \\ j &= 2, \dots, n \\ & \quad u_{1j} = \frac{a_{1j}}{l_{11}} \\ j &= 2, \dots, n-1 \\ & \quad l_{jj} = a_{jj} - \sum_{k=1}^{j-1} l_{jk} u_{kj} \\ & \quad i = j+1, m \\ & \quad \quad l_{ij} = a_{ij} - \sum_{k=1}^{j-1} l_{ik} u_{kj} \\ & \quad \quad u_{ji} = \frac{a_{ji} - \sum_{k=1}^{j-1} l_{jk} u_{ki}}{l_{jj}} \\ l_{nn} &= a_{nn} - \sum_{k=1}^{n-1} l_{nk} u_{kn} \end{align*} $$

Seudocódigo

function lu_crout(A)
    for i=1 to m do
        l(i,1) = a(i,1)
    end for
    for j=2 to n do
        u(1,j) = a(1,j)/l(1,1)
    end for
    for j=2 to n-1 do
        suma = a(j,j)
        for k=1 to j-1 do
            suma = suma - l(j,k)*u(k,j)
        end for
        l(j,j) = suma
        for i=j+1 to m do
            sumav = a(i,j)
            sumah = a(j,i)
            for k=1 to j-1 do
                sumav = sumav - l(i,k)*u(k,j)
                sumah = sumah - l(j,k)*u(k,i)
            end for
            l(i,j) = sumav
            u(j,i) = sumah/l(j,j)
        end for
    end for
    suma = a(n,n)
    for k=1 to n-1 do
        suma = suma - l(n,k)*u(k,n)
    end for
    l(n,n) = suma
    return a
end function

Implementación descomposición de Crout

import numpy as np

def lu_crout(a):
    m, n = a.shape
    u = np.eye(m)
    l = np.zeros((m,n))
    for i in range(0,m):
        l[i,0] = a[i,0]

    for j in range(1,n):
        u[0,j] = a[0,j]/l[0,0]

    for j in range(1,n-1):
        suma = a[j,j]
        for k in range(0,j):
            suma = suma - l[j,k]*u[k,j]
        l[j,j] = suma
        for i in range(j+1,m):
            sumav = a[i,j]
            sumah = a[j,i]
            for k in range(0,j):
                sumav = sumav - l[i,k]*u[k,j]
                sumah = sumah - l[j,k]*u[k,i]
            l[i,j] = sumav
            u[j,i] = sumah/l[j,j]

    suma = a[n-1,n-1]
    for k in range(0,n-1):
        suma = suma - l[n-1,k]*u[k,n-1]
    l[n-1,n-1] = suma

    return l, u

A = np.array([[1,1,2,3],
              [2,1,-1,1],
              [3,-1,-1,2],
              [-1,2,3,-1]],float)
print(A)

    [[ 1.  1.  2.  3.]
     [ 2.  1. -1.  1.]
     [ 3. -1. -1.  2.]
     [-1.  2.  3. -1.]]

#resultado
L, U = lu_crout(A)
print(L)
print(U)

    [[  1.   0.   0.   0.]
     [  2.  -1.   0.   0.]
     [  3.  -4.  13.   0.]
     [ -1.   3. -10.  -3.]]
    [[ 1.  1.  2.  3.]
     [ 0.  1.  5.  5.]
     [ 0.  0.  1.  1.]
     [ 0.  0.  0.  1.]]

#revisando el resultado
np.dot(L,U)

    array([[ 1.,  1.,  2.,  3.],
           [ 2.,  1., -1.,  1.],
           [ 3., -1., -1.,  2.],
           [-1.,  2.,  3., -1.]])

A = np.array([[3,-0.1,-0.2],
              [0.1,7,-0.3],
              [0.3,-0.2,10]],float)
print(A)

    [[  3.   -0.1  -0.2]
     [  0.1   7.   -0.3]
     [  0.3  -0.2  10. ]]

#resultado
L, U = lu_crout(A)
print(L)
print(U)

    [[  3.           0.           0.        ]
     [  0.1          7.00333333   0.        ]
     [  0.3         -0.19        10.01204188]]
    [[ 1.         -0.03333333 -0.06666667]
     [ 0.          1.         -0.04188482]
     [ 0.          0.          1.        ]]

#revisando el resultado
np.dot(L,U)

    array([[  3. ,  -0.1,  -0.2],
           [  0.1,   7. ,  -0.3],
           [  0.3,  -0.2,  10. ]])