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

4.2 Resolución de varios sistemas lineales

En ocasiones podemos tener varios sistemas lineales en los que la matriz del sistema [math]A[/math] es siempre la misma, y en los que el vector de términos independiente [math]b_k \ \ \ k\in[1,M][/math] es diferente para cada sistema.

Los métodos exactos para resolver sistemas numéricamente realizan del orden de [math]N^3[/math] operaciones para resolver un sistema, donde [math]N[/math] es el orden de la matriz [math]A[/math]. Por tanto, para resolver [math]M[/math] sistemas necesitaríamos [math]MN^3[/math] operaciones en total. El método de Doolittle también necesita [math]N^3[/math] operaciones para obtener los factores. Pero la resolución de sistemas lineales solo necesita del orden de [math]N^2[/math] operaciones.

Por tanto, si usamos el método de Doolittle para obtener los factores de la matriz [math]A[/math], una vez hallados, podremos resolver cada uno de los sistemas triangulares siguientes

1. [math]Ly_k = b_k[/math]
2. [math]U x_k = y_k[/math]

En total se necesitarán [math]N^3 + 2MN^2[/math] operaciones para resolver los [math]M[/math] sistemas. Si [math]M[/math] es grande, esta opción conlleva muchas menos operaciones que [math]MN^3[/math].

4.3 Cálculo de la matriz inversa

El cálculo de la matriz inversa es similar al caso anterior de resolver varios sistemas lineales con la misma matriz del sistema [math]A[/math]. En este caso, los vectores de términos independientes son cada una de las columnas de la matriz inversa [math]i_k[/math]. Cada una de las soluciones [math]x_k[/math] obtenidas de esta manera es la columna correspondiente de la matriz inversa de [math]A[/math].

5 Ventajas del método de Doolittle

Los ordenadores cuenta con diferentes elementos como procesador, memoria, disco duro, puertos USB, unidades de DVD. En el cálculo numérico, los dos elementos fundamentales que deciden si el cálculo se realizará antes o después son el procesador y la memoria.

El procesador se encarga de realizar las diferentes operaciones matemáticas, y la memoria contiene los valores de las variables y de los resultados de las operaciones. Existen diferentes maneras de organizar el procesador y la memoria en un ordenador. A las diferentes maneras de organizar (y al tipo de operaciones que puede realizar el procesador) se las conoce como arquitectura.

Arquitectura de von Neumann

Los ordenadores modernos se organizan todos usando la denominada arquitectura de von Neumann. En esta arquitectura, la memoria y el procesador están conectados por medio de un bus, que se encarga de pasar la información entre el ordenador y la memoria. El bus actúa de cuello de botella en el sistema. Por muy rápido que sea el procesador, si no puede recibir los datos que necesita para trabajar al ritmo adecuado, entonces la velocidad de cómputo estará limitada.

Para solucionar los problemas del cuello de botella del bus en los ordenadores más recientes existen dos tipos de memoria: la memoria principal y la caché. La caché está conectada al procesador por un bus mucho más ancho y rápido. A su vez, la caché organiza la información para ir requiriendo datos de la memoria principal y mantener siempre datos en la caché para que el procesador los pueda usar sin tener que esperar a que estén disponibles.