viernes, 30 de septiembre de 2016

Factorización LU

En esta ocasión estudiamos la factorización $LU$ y la eficiencia que promueve en la solución de sistemas de ecuaciones con muchas incógnitas. Como caso particular, se aplica en el método de diferencias finitas que tiene lugar cuando se quiere resolver algunas ecuaciones diferenciales lineales de orden 2 con valores en la frontera. Esperamos que sea de su agrado.

Una de las factorizaciones más elegante e importante del Álgebra Lineal es la factorización $A = LU$ Cuando se habla de la posibilidad de resolver un sistema de ecuaciones por un método más eficiente que la eliminación Gaussiana cuesta creerlo! Pero existe dicho método y esto es gracias a la factorización de la matriz $A$ como producto de dos matrices $L$ (triangular inferior) y $U$ (triangular superior).

Para entender la factorización $LU$ de una matriz $A$ a la que se puede aplicar eliminación Gaussiana sin encontrar ceros en la diagonal principal y sólo con operaciones elementales de fila que consisten en sustituir una fila, por el resultado de sumarle a ésta un múltiplo escalar de otra, una estrategia es pensar en dicha operación elemental como una premultiplicación por una matriz, en este caso, una matriz elemental triangular inferior.

Suponga que durante el proceso de eliminación Gaussiana para la matriz $$A = \left [ \begin{array}{ccc} 2 & 2& 2\\ 4& 7& 7\\6 &18 &22 \end{array} \right ] $$ A las operaciones elementales de fila:

  1. Sustituir fila $2$ por fila $2$ $-$ $2$ veces fila $1$.
  2. Sustituir fila $3$ por fila $3$ $-$ $3$ veces fila $1$.

Se les interpreta como un premultiplicación por la matriz $T = \left [ \begin{array}{ccc} 1 &0 & 0\\-2 & 1 & 0\\ -3& 0 & 1 \end{array} \right ] $. El cálculo directo de $T$ por $A$ es: $$\left [ \begin{array}{ccc} 1 &0 & 0\\-2 & 1 & 0\\ -3& 0 & 1 \end{array} \right ] \left [ \begin{array}{ccc} 2 & 2& 2\\ 4& 7& 7\\6 &18 &22 \end{array} \right ] = \left [ \begin{array}{ccc} 2 & 2& 2\\ 0& 3& 3\\0 &12 &16 \end{array} \right ] $$ El resultado obtenido es el mismo que si se aplican las operaciones elementales de fila.
Pero, cómo se construye en general dicha matriz $T$?, además, es $T$ no singular (invertible)?. Esas preguntas se van a responder y se vará la relación que tienen con la factorización $LU$.

Para ver que lo anterior (en el ejemplo ilustrativo inicial) sigue pasando en general, se define una matriz elemental triangular interior de orden $n$ como $$T _{k} = I - c _{k} e _{k} ^{T} $$ donde $c _{k} $ es columna con $k$ ceros al principio, es decir de la forma $c _{k} = \left [ \begin{array}{c} O _{k} \\ -\mu _{k + 1 } \\ \vdots \\ -\mu _{n} \end{array} \right ] $ y $e _{k} = I _{*k} $ (k-ésima columna de la identidad $I _{n \times n}$).
La forma que tiene $T _{k} $ es similar a la identidad $I$, excepto que debajo del 1 en la $k$-ésima columna tiene los valores $-\mu _{k + 1 },..., -\mu _{n} $, de la siguiente forma: $$T _{k} = \left [ \begin{array}{cccccc} 1 _{1} & & & & & \\ & \ddots & & & & \\ & & 1 _{k} & & & \\ & &-\mu _{k+1} &1 _{k+1} & &\\& &\vdots & &\ddots &\\& & -\mu _{n} & & & 1 _{n} \end{array} \right ] $$ Considerando los espacios vacíos como ceros, note que tiene la misma forma que la matriz $T$ en el ejemplo ilustrativo inicial.

Se puede sospechar que $T$ sea invertible, pero lo que ahora es sorprendente es que $T _{k} ^{-1} = I + c _{k} e _{k} ^{T}$, lo que se puede ver por multiplicación directa $$(I + c _{k} e _{k} ^{T} )(I - c _{k} e _{k} ^{T} ) = I - c _{k} e _{k} ^{T} + c _{k} e _{k} ^{T} - c _{k} e _{k} ^{T} c _{k} e _{k} ^{T} $$ No queda claro aún es por qué el término $c _{k} e _{k} ^{T} c _{k} e _{k} ^{T} $ es la matriz nula. Para ello, note que $$c _{k} e _{k} ^{T} c _{k} e _{k} ^{T} = c _{k} ( e _{k} ^{T} c _{k} )e _{k} ^{T} = c _{k} 0 e _{k} ^{T} $$ Se ha sustituido $e _{k} ^{T} c _{k} $ por el n\'umero $0$ puesto que recuerde que $e _{k} ^{T} $ solo tiene un 1 en la posición $k$, mientras que $c _{k}$ tiene un $0$ en ese mimo lugar.
Se puede decir que $T _{k} ^{- 1} $ también es una matriz elemental triangular inferior. $T _{k} ^{-1} $ es muy similar a $T _{k} $, los números fuera de la diagonal de 1's cambian de signo.

Todo lo anterior da una idea de qué es lo que se hace al poner ceros bajo la diagonal en el proceso de escalonar $A$ para llegar a una matriz $U$ triangular superior. En realidad lo que se hace es premultiplicar por $n - 1$ matrices elementales triangulares inferiores, es decir: $$A \to T _{1} A \to T _{2} T _{1} A \to \cdots \to T_{n - 1 } T _{n - 2} T _{n - 3} \cdots T _{2} T _{1} A = U $$ Y al premultiplicar por las inversas $$T _{n - 2}\cdots T _{2} T _{1} A = T _{n - 1 } ^{-1} U \ \to \ T _{n - 3}\cdots T _{2} T _{1} A = T _{n - 2} ^{-1} T _{n - 1 } ^{-1} U \ \to \ A = T _{1} ^{-1} \cdots T _{n - 1} ^{-1} U $$ Asi que $$A = LU$$ con $L = T _{1} ^{-1} \cdots T _{n - 1} ^{-1} $.

Ahora, $L =T _{1} ^{-1} \cdots T _{n - 1} ^{-1} = (I + c _{1} e _{1} ^{T} )(I + c _{2} e _{2} ^{T} )(I + c _{3} e _{3} ^{T} ) \cdots (I + c _{n-1} e _{n-1} ^{T} ) $. Por multiplicación directa se puede notar que $$L = I + c _{1} e _{1} ^{T} + c _{2} e _{2} ^{T} + c _{3} e _{3} ^{T} + \cdots+ c _{n-1} e _{n-1} ^{T} $$ Puesto que en el producto de las matrices del tipo $I + c _{k} e _{k} ^{T} $ se cancelan muchos térmidos dado que los factores del tipo $ e _{k} ^{T} c _{j} $ son el número cero cuando $k < j$. (Recuerde que $e _{k} ^{T}$ solo tiene un 1 en la posición $k$ mientras que $c _{j} $ tiene posiblemente números diferentes de cero sólo hasta la posición $j + 1$).

Para ver cómo es $L$ se debe percatar de la forma de la suma $$I + c _{1} e _{1} ^{T} + c _{2} e _{2} ^{T}+ c _{3} e _{3} ^{T} + \cdots+ c _{n-1} e _{n-1} ^{T} $$ $$I +\left [ \begin{array}{cccccc} & & & & & \\ l _{21} & & & & & \\ l _{31}& & & & & \\ l _{41} & && & &\\ \vdots& & & & &\\ l _{n1} & & & & & \end{array} \right ]+ \left [ \begin{array}{cccccc} & & & & & \\ & & & & & \\&l _{32} & & & & \\ & l _{42} && & &\\ & \vdots & & & &\\ & l _{n2}& & & & \end{array} \right ] + \cdots + \left [ \begin{array}{cccccc} & & & & & \\ & & & & & \\& & & & & \\ & && & &\\ & & & & &\\ & & & & l _{n,n - 1} & \end{array} \right ] $$ Es decir $$L = \left [ \begin{array}{cccccc} 1 & & & & & \\ l _{21} & 1& & & & \\ l _{31}&l _{32} & 1 & & & \\ l _{41} & l _{42} &l _{43} & \ddots & &\\\vdots & \vdots & \vdots & &1 &\\ l _{n1} & l _{n2}& l _{n3}& & l _{n,n - 1} & 1 \end{array} \right ] $$

A modo de conclusión se tiene que:

  1. $L$ es triangular inferior y $U$ es triangular superior. Además, $l _{ii} = 1 $ y $u _{ii} \neq 0 $ para cada $i = 1,2,..., n$.
  2. Bajo la diagonal de $L$, la entrada $l _{ij} $ es el múltiplo de la fila $j$ que es sustraida de la fila $i$ con el objetivo de anular la posición $ij$ durante la eliminación Gaussiana.
  3. $U$ es el resultado final al escalonar $A$. Además, las matrices $L$ y $U$ son únicas.

En el ejemplo inicial con $$A = \left [ \begin{array}{ccc} 2 & 2& 2\\ 4& 7& 7\\6 &18 &22 \end{array} \right ] $$ se puede terminar la eliminación y notar que $$L= \left [ \begin{array}{ccc} 1 & 0& 0\\ 2& 1& 0\\3 &4 &1 \end{array} \right ] \ \ \ y \ \ \ A = \left [ \begin{array}{ccc} 2 & 2& 2\\ 0& 3& 3\\0 &0 &4 \end{array} \right ] $$ Es importante percatarse de algo genial, y es que las entradas que caracterizan a $L$ (sólo las de la parte inferior izquierda) y las que caracterizan a $U$ (sólo las de la parte superior derecha), suman tantas entradas como las que inicialmente tiene $A$, lo que significa que se necesita la misma capacidad de memoria para almacenar sólo a $A$ que para almacenar a las dos matrices $L$ y $U$. Ahora se da paso al análisis del proceso con el que se resuelve un sistema de ecuaciones $Ax = b$ usando $L$ y $U$ en lugar de $A$. Además, se resalta la importancia de tener las componentes $L$ y $U$ para resolver familias de sistemas $Ax = \overline{b}$ que tengan la misma matriz de coeficientes $A$ y cualquier otro vector de términos independientes $\overline{b}$.
Suponga que ya se tienen $L$ y $U$. El sistema $Ax = b$ se puede ver como $LUx = b$ y si se asume $Ux = y$ entonces el sistema inicial $Ax = b$ se puede interpretar como dos sistemas triangulares $Ly = b$ y $Ux = y $.
Se procede con el sistema $Ly = b$ en donde el vector $y$ de incógnitas se halla por sustitución directa $$ \left [ \begin{array}{cccccc} 1 & & & & & \\ l _{21} & 1& & & & \\ l _{31}&l _{32} & 1 & & & \\ l _{41} & l _{42} &l _{43} & \ddots & &\\\vdots & \vdots & \vdots & &1 &\\ l _{n1} & l _{n2}& l _{n3}& & l _{n,n - 1} & 1 \end{array} \right ] \left [ \begin{array}{c} y _{1} \\ y _{2} \\ y _{3} \\ y _{4}\\ \vdots \\ y _{n} \end{array} \right ] = \left [ \begin{array}{c} b _{1} \\ b _{2} \\ b _{3} \\ b _{4} \\ \vdots \\ b _{n} \end{array} \right ] $$ Es decir, $y _{1} = b _{1}$, $\ \ y _{2} = b _{2} - l _{21} y _{1}$, $\ \ y _{3} = b _{3} - l _{31} y _{1} $, etc.
En otras palabras, $y _{1} = b _{1} \ \ $ y $\ \ y _{i} = b _{i} - \displaystyle \sum_{k = 1 } ^{i - 1 } l _{ik} y _{k} $ para $i = 2,3,...,n$.
Cuando ya se tienen las entradas del vector $y$, se monta el sistema triangular superior $Ux = y$ que se resuelve usando sustitución regresiva comenzando con $x _{n} = \frac{y _{n} }{u _{n} } $ y tomando $$ x _{i} = \frac{1 }{u _{ii} } \left ( y _{i} - \displaystyle \sum_{k = i + 1 } ^{n} u _{ik} x _{k} \right ) \ \ \ para \ \ \ i = n - 1 , n - 2, ..., 1.$$ Cabe notar que con éste método de sustituciones usando $L$ y $U$, sólo se requieren $n ^{2} $ multiplicación y $n ^{2} - n$ adiciones. Por otro lado, para la eliminación en el sistema aumentado $A | b$ se requieren alrededor de $n ^{3} /3$ multiplicaciones, lo que indica que para $n$ grande, usar la factorización $LU$ genera evidente eficiencia. Además, en cualquier otro sistema $Ax = \overline{b}$ con la misma matriz $A$ se pueden aplicar los mismo algoritmos eficientes de sustitución.
Por ejemplo, en el importante problema (atacado con método de diferencias finitas) para la ecuación diferencial $\frac{d ^{2} }{dt ^{2} } x = f(t) ,\ \ \ con \ \ \ 0 \leq t \leq 1, \ \ x(0) = x(1) = 0$, sale a flote la matriz de coeficientes $$A = \left [ \begin{array}{ccccccc} 2& -1&& && & \\ -1& 2&-1& && & \\ &-1 &2& -1&& &\\ & &&\ddots && &-1 \\ & && &&-1 &2 \end{array} \right ] $$ Cuyos factores son $$L = \left [ \begin{array}{ccccccc} 1& && && & \\ \frac{-1}{2} & 1& && & & \\ &\frac{-2}{3} &1& && &\\ & &\frac{-3}{4}&\ddots && & \\ & && && \frac{-(n-1)}{n}&1 \end{array} \right ] \ \ \ y \ \ \ U = \left [ \begin{array}{ccccccc} 2& -1&& && & \\ &\frac{3}{2} &-1& && & \\ & &\frac{4}{3}& -1&& &\\ & &&\ddots && &-1 \\ & && && &\frac{n + 1 }{n} \end{array} \right ] $$ Y cuyas fórmulas de sustitución recursiva para los respectivos sistemas triangulares $Ly = b$ y $Ux = y$ son respectivamente $$y _{1} = b _{1} , y _{k} = b _{k} + \frac{k - 1 }{k} y _{k - 1 }, \ para \ k = 2,...,n. $$ $$x _{n} = \frac{n}{n + 1 }y _{n}, \ x _{k} = \frac{k}{k + 1 } (y _{k} + x _{k+1} ),\ para \ \ k = n - 1, n - 2,...,1. $$

Se recomienda ir a los trabajos de Carl D Meyer y de Gilbert Strang

Referencias

0 comentarios: