Factorización de Doolittle

Este artículo está en versión beta. El autor de este artículo no lo ha terminado todavía, por favor no lo edites hasta que elimine este mensaje.

En álgebra lineal, se conoce por factorización de matrices al proceso que a partir de una matriz cuadrada [math]A[/math] halla dos matrices triangulares inferior y superior, tal que [math]A = L U[/math], donde [math]L[/math] es la matriz triangular inferior (L de lower) y [math]U[/math] es la matriz triangular superior (U de upper). Existen muchos métodos numéricos para obtener estas matrices [math]L[/math] y [math]U[/math], y su obtención tiene aplicaciones en la resolución de sistemas lineales, cálculo de determinantes y en el cálculo de matrices inversas. En este artículo estudiamos el método de Doolittle para obtener estas matrices. Este método tiene la particularidad que hace que la diagonal de la matriz [math]L[/math] sea unitaria.

1 Métodos de factorización

Para resolver sistemas de ecuaciones lineales numéricamente existen dos grandes familias de métodos:

• Métodos exactos
• Métodos iterativos

Los métodos exactos proporcionan una solución exacta del sistema. A pesar de ser considerados métodos numéricos, el procedimiento no es una aproximación a la solución, sino la solución en sí. En la mayoría de las ocasiones es preferible usar estos métodos. Sin embargo, en ocasiones no es posible aplicarlos (por ejemplo, con matrices mal condicionadas), o puede ser más costoso usarlos por las propiedades de las matrices (por ejemplo, matrices sparse, que contienen muchos ceros). En estos casos, un método iterativo proporciona una solución aproximada del sistema.

En cuanto a los métodos exactos para factorizar una matriz, existen también diferentes tipos de métodos. El método de Doolittle es un método denominado compacto, porque es sencillo de programar y requiere poca memoria una vez implementado en el ordenador. Existen muchos métodos compactos de factorización de matrices, y suelen diferir en el tratamiento que hacen de los elementos de la diagonal de las matrices [math]L[/math] y/o [math]U[/math]. En concreto, el método de Doolittle genera una matriz [math]L[/math] que tiene [math]1[/math] en todos los elementos de la diagonal. Este método proporciona una ventaja de cálculo en ordenadores modernos, que cuentan con memoria caché además de la memoria principal. En este tipo de arquitecturas de ordenador, la secuencia en la que se realizan los cálculos puede ser más importante que la cantidad de cálculos que se realizan. El método de Doolittle tiene una secuencia de operaciones óptima para ejecutarse en un ordenador con memoria caché.

2 Descripción del método

El resultado del método serán dos matrices triangulares tales que [math]A=LU[/math]. La matriz [math]L[/math] será triangular inferior y de diagonal unitaria. La matriz [math]U[/math] será triangular superior. Al ser las matrices triangulares, podemos calcular cada elemento de la matriz [math]A[/math] mediante

[math]\displaystyle a_{ij} = \sum_{k=1}^{\min(i,j)} l_{ik} u_{kj}[/math]

donde [math]a[/math], [math]l[/math] y [math]u[/math] se refieren a los elementos de las matrices [math]A[/math], [math]L[/math] y [math]U[/math], respectivamente. El límite superior de la suma se debe a que las matrices triangulares tienen ceros en su mitad inferior o superior, y por tanto, el producto matricial sumaría cero.

En el caso de que [math]i \leq j[/math], es decir, el número de fila es más pequeño que el número de la columna, en otras palabras, estamos por encima de la diagonal (o en la diagonal), entonces [math]\min(i,j)= i[/math] y tendremos

[math]\displaystyle l_{ii}u_{ij} = u_{ij} = a_{ij} - \sum_{k=1}^{i-1} l_{ik} u_{kj}[/math]

ya que [math]l_{ii} = 1[/math]. Es decir, hemos despejado del sumatorio el sumando cuando [math]k=i[/math], y eso nos permite obtener una expresión genérica para calcular [math]u_{ij}[/math].

Si empezamos en [math]i=1[/math] la ecuación anterior nos dice que [math]u_{1j} = a_{1j}[/math], por lo que ya hemos calculado la primera fila de la matriz [math]U[/math].

Para seguir calculando filas de [math]U[/math] es necesario conocer los valores de [math]L[/math]. Si nos movemos por debajo de la diagonal, [math]j\leq i[/math], tendremos que [math]\min(i,j)=j[/math], y obtenemos

[math]\displaystyle l_{ij} = \frac{a_{ij} - \sum_{k=1}^{j-1} l_{ik} u_{kj}}{u_{jj}}[/math]

Por tanto, podemos calcular [math]\displaystyle l_{i1} = \frac{a_i}{u_{11}}[/math], es decir, hemos calculado ya la primera columna de [math]L[/math]. Al conocer la primera columna de [math]L[/math] podemos intentar calcular la siguiente fila de [math]U[/math]. Al terminar una fila de [math]U[/math], ya conocemos los elementos para calcular la siguiente columna de [math]L[/math]. Aplicando sucesivamente este método, obtendremos de manera completa las matrices [math]L[/math] y [math]U[/math].

3 Implementación en código MATLAB

Debido a que las ecuaciones del método recorren matrices en orden de filas y columnas, es relativamente sencillo implementar este método en código MATLAB. Veamos el código de una función que devuelve las dos matrices [math]L[/math] y [math]U[/math] a partir de una matriz cuadrada [math]A[/math]:

function [L,U] = doolittle (A) 
% Reservar espacio para las matrices L y U
% Empezando con ceros hay sitios que ya no tendremos que tocar
L = zeros(size(A));
U = zeros(size(A));
[nRows, nColumns] = size(A);
% Recorremos en orden de columnas, es más rápido en MATLAB
for i=1:nRows
  for j=1:nColumns    
    % Estamos por encima de la diagonal, hallamos elemento de U
    if i<=j
      U(i,j) = A(i,j);
      for k=1:i-1
        U(i,j) = U(i,j) - L(i,k)*U(k,j);
      end
    end    
    % Estamos por debajo de la diagonal, hallamos elemento de L
    if j<=i
      L(i,j) = A(i,j);
      for k=1:j-1
        L(i,j) = L(i,j) - L(i,k)*U(k,j);
      end
      L(i,j) = L(i,j)/U(j,j);
    end
  end
end

4 Aplicaciones del método

4.1 Resolución de un sistema lineal

La aplicación más directa de la factorización LU es la resolución de sistemas lineales. Al transformar la matriz del sistema en un producto de matrices triangulares, el sistema de ecuaciones puede transformarse en dos sistemas triangulares, que son muy sencillos de resolver en un ordenador.

Si tenemos el sistema de ecuaciones [math]Ax = b[/math], podemos cambiarlo por [math]LUx=b[/math]. Pero [math]Ux[/math] será también un vector columna, por lo que nos queda por un lado que [math]Ly=b[/math], y por otro que [math]Ux=y[/math]. Ambos son sistemas triangulares. Resolvemos el primero para obtener [math]y[/math], y a continuación el segundo para obtener la solución buscada [math]x[/math].