Factorización de Doolittle

De MateWiki
Saltar a: navegación, buscar

En álgebra lineal, se conoce por factorización LU de matrices[1] 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 resolución de sistemas lineales

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, pueden intentar resolver el sistema directamente, o usando sus factores. El método de Doolittle es un método de factorización LU 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é.

Los métodos de factorización LU se pueden usar para resolver sistemas lineales, entre otras aplicaciones.

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{\displaystyle 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. La implementación típica requiere recorrer una matriz, en este caso recorremos siempre empezando por columnas, para mejorar el rendimiento. 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]:

 1 function [L,U] = doolittle (A) 
 2 % Reservar espacio para las matrices L y U
 3 % Empezando con ceros hay sitios que ya no tendremos que tocar
 4 L = zeros(size(A));
 5 U = zeros(size(A));
 6 [nRows, nColumns] = size(A);
 7 % Recorremos en orden de columnas, es más rápido en MATLAB
 8 for j=1:nColumns
 9   for i=1:nRows
10     % Estamos por encima de la diagonal, hallamos elemento de U
11     if i<=j
12       U(i,j) = A(i,j);
13       for k=1:i-1
14         U(i,j) = U(i,j) - L(i,k)*U(k,j);
15       end
16     end    
17     % Estamos por debajo de la diagonal, hallamos elemento de L
18     if j<=i
19       L(i,j) = A(i,j);
20       for k=1:j-1
21         L(i,j) = L(i,j) - L(i,k)*U(k,j);
22       end
23       L(i,j) = L(i,j)/U(j,j);
24     end
25   end
26 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 usando el algoritmo de sustitución hacia atrás.

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

Arquitectura de von Neumann (el bus es un cuello de botella)
Arquitectura de von Neumann con jerarquía de memorias (caché y memoria principal, aliviando el cuello de botella)

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.

Los ordenadores modernos se organizan todos usando la denominada arquitectura de von Neumann[2]. 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. En esta arquitectura, toda la información que se transfiere entre los diferentes elementos del ordenador pasa a través del bus, y se comunica a todos los elementos. Cada elemento toma la información que va dirigida a él, e ignora cualquier información dirigida a otro elemento.

Hace décadas, dada la velocidad de proceso de los ordenadores, este sistema era una buena idea, ya que el bus nunca dejaba ocioso el procesador del ordenador. Sin embargo, en los ordenadores actuales, el bus actúa de cuello de botella en el sistema, ya que el procesador consume información a mayor ritmo del que el bus es capaz de suministrarla. 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. Es imposible aprovechar toda la velocidad de cómputo del procesador.

Para solucionar los problemas del cuello de botella del bus[3] 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. Esta jerarquía de memorias mejora el rendimiento, sobre todo en aplicaciones gráficas de escritorio. No se pensó como una mejora para potenciar el cálculo numérico, pero puede aprovecharse para elegir métodos optimizados para este tipo de arquitecturas.

Tradicionalmente, la diferencia entre la velocidad de cómputo con un método u otro radicaba en la complejidad del método, es decir, en la cantidad de operaciones que tiene que realizar. Sin embargo, al introducir una memoria caché, que tiene una comunicación mucho más rápida con el procesador, la secuencia de operaciones es ahora también un parámetro fundamental para determinar la velocidad de cómputo.

Al ser la memoria caché normalmente más pequeña que la memoria principal (va instalada en la misma placa que el procesador, por lo que su tamaño es muy limitado), no se puede guardar toda la información en la memoria caché. Por ejemplo, típicamente estas memorias son de 1MB o 2MB, por lo que es imposible almacenar las matrices [math]A[/math], [math]L[/math] y [math]U[/math] en esa memoria siempre. Por ello, nuestro método de cálculo tiene que estar pasando información entre la caché y la memoria principal.

En estas situaciones, si el procesador tiene siempre en la caché los valores que necesita, irá más rápido que si tiene que ir pidiendo valores a la memoria principal. Por la manera en que funciona esta arquitectura, la caché siempre contendrá los valores necesarios para una operación, y los valores resultantes de la operación. Esta característica es precisamente la que explota el método de Doolittle. Cada vez que calculamos un valor de [math]u_{ij}[/math], podemos calcular el valor de [math]l_{ij}[/math] correspondiente. Pero para calcular el siguiente valor de [math]u_{ij}[/math] necesitamos el valor de [math]l_{ij}[/math] que acabamos de calcular. Gracias a esta característica del método, se minimizan las peticiones de datos a la memoria principal, y aunque el número de operaciones es el mismo, su ejecución suele ser más rápida. Las ganancias de tiempo de ejecución son del orden del 50% en el mejor de los casos.

Esta ventaja es más difícilmente explotable desde MATLAB u Octave UPM, ya que son lenguajes de muy alto nivel, y es imposible controlar si los datos se guardan en caché o no. Además, al ser entornos de programación, realizan muchas otras operaciones ocultas que también requieren del uso de la caché, por lo que no está claro que se aproveche totalmente la caché con este método. En cualquier caso, por la manera en la que funcionan MATLAB y Octave podemos facilitar la reutilización de la caché recorriendo siempre las matrices primero en columnas, tal y como se hace en la implementación del método mostrada arriba.

6 Referencias

  1. LU decomposition (Wikipedia EN)
  2. von Neumann architecture (Wikpiedia EN)
  3. von Neumann bottleneck (Wikipedia EN)