Factorización gaussiana LU de una matriz regular

Preliminares

Toda matriz de orden $n$, regular, $A$, cuyos determinantes principales sean distintos de cero, puede expresarse de la forma $A=LU$ mediante operaciones de reducción gaussiana sin necesidad de pivotaje, donde $L$ es una matriz triangular inferior de orden $n$ con unos en la diagonal principal y $U$ es una matriz triangular superior de orden $n$, siendo ambas matrices regulares (sus determinantes son no nulos).

En estas condiciones, las matrices $L$ y $U$ son únicas. En efecto, supongamos que no lo sean, esto es, si podemos escribir $A=L_1\,U_1=L_2\,U_2$; entonces, como todas estas matrices son inversibles, podemos escribir $L_{2}^{-1}\,L_1\,U_1=L_{2}^{-1}\,L_2\,U_2$, luego $L_{2}^{-1}\,L_1\,U_1=I\,U_2$ y por tanto $L_{2}^{-1}\,L_1\,U_1=U_2$; multiplicando ahora por $U_{2}^{-1}$ por la derecha en ambos miembros, tenemos $L_{2}^{-1}\,L_1\,U_1\,U_{2}^{-1}=U_2\,U_{1}^{-1}$, esto es $L_{2}^{-1}\,L_1\,U_1\,U_{1}^{-1}=U_2\,U_{1}^{-1}$ y por tanto, $L_{2}^{-1}\,L_1\,I=U_2\,U_{1}^{-1}$, que es lo mismo que $L_{2}^{-1}\,L_1=U_2\,U_{1}^{-1}$. La matriz inversa de una matriz triangular inferior es otra matriz triangular inferior, y lo mismo sucede con la inversa de una matriz triangular superior, así que, como el producto de dos matrices triangulares inferiores es otra matriz triangular inferior (primer miembro) y el producto de dos matrices triangulares superiores es una matriz triangular superior, la única manera de que esto sea posible es que ambos miembros de la igualdad sean la matriz identidad, por consiguiente: $L_{2}^{-1}\,L_1=I \Leftrightarrow L_2=L_1$ y $U_2\,U_{1}^{-1}\Leftrightarrow U_2=U_1$. $\square$

Observaciones:

• En una futura entrada trataré el caso en que sea necesario el pivotaje por el hecho de encontrar elementos nulos en la diagonal de alguno de los pasos de reducción y que supondrían divisiones por cero —sería el caso en el que no todos los determiantes principales sea distinto de cero, incumpliendo la condición dada arriba—, o bien números muy pequeños en algún elemento de la diagonal que harían que los errores de redondeo se amplificasen enormemente. Trataremos pues más adelante los casos de pivotaje, bien sea pivotaje parcial (también llamado maximal por columnas) en el que hay que permutar filas de la matriz; o bien se trate de pivotaje total (también llamado pivotaje completo), en el que no sólo hay que permutar filas sino también columnas. En los casos de pivotaje veremos que intervendran las llamadas matrices de permutación.
• Si no son necesarias operaciones de pivotaje, hay otros métodos para expresar una matriz $A$ como producto de una matriz triangular inferior por una matriz triangular superior, pero que no se basan en la reducción gaussiana:
  • El método de Doolittle permite una factorización de este tipo, obteniendo la matriz triangular inferior $L$ con unos en la diagonal principal, sin utilizar algoritmos de reducción (véase este ejemplo práctico).
  • El método de Crout permite también una factorización similar, obteniendo la matriz triangular superior $U$ con unos en la diagonal principal, sin utilizar algoritmos de reducción (véase este ejemplo práctico)
• Además de la factorización de $A$ mediante el producto de una matriz triangular inferior por una matriz triangular superior hay otras formas de factorización:
  • Si $A$ es una matriz simétrica, puede realizarse una factorización más conveniente de la forma $A=LDL^\top$, donde $L$ es una matriz triangular inferior con unos en la diagonal y $D$ es una matriz diagonal.
  • Si $A$ es una matriz definida positiva, el criterio de Sylvester asegura que los elementos de la diagonal, $d_{kk}$, son positivos. Y definiendo $\mathcal{L}:=L\,D^{\frac{1}{2}}$, con $D^{\frac{1}{2}}=\text{diagonal}(d_{11}^{\frac{1}{2}},\ldots,d_{nn}^{\frac{1}{2}})$, entonces la matriz $A$ puede expresarse de la forma $A=\mathcal{L}\,\mathcal{L}^\top$, que se denomina factorización de Cholesky.
    • Si $A$ es simétrica y definida positiva, entonces la factorización de Cholesky de $A$ que se acaba de exponer es más sencilla, pues $\mathcal{L}=L$. Así, $A=L\,L^\top$; en otras palabras, en tales condiciones la matriz triangular superior $U$ de la factorización $LU$ es $U=L^\top$

-oOo-
Ensayemos ahora un poco con una matriz regular como, por ejemplo, $A=\begin{pmatrix}3&-1&1\\1&1&1\\2&1&0\end{pmatrix}$, cuya reducción no requiere pivotamiento parcial pues no nos vamos a encontrar con ceros en la diagonal principal de la matrices que iremos obteniendo.

Más adelante, extenderemos lo que vamos a deducir ahora al caso en que sí se requiere pivotamiento.

Empezaremos pues con el proceso de reducción. Siendo esta matriz de orden $3$, tendremos que llevar a cabo dos etapas.

Obtenemos ceros por debajo de la columna del elemento $a_{11}$ haciendo las siguientes transformaciones $\begin{matrix}f_1 \rightarrow f_1 \\ -\frac{1}{3}f_1+f_2\rightarrow f_2\\-\frac{2}{3}f_1+f_3\rightarrow f_3\end{matrix}$, obteniendo la matriz $\begin{pmatrix}3&-1&1\\0&\frac{4}{3}&\frac{2}{3}\\0&\frac{5}{3}&-\frac{2}{3}\end{pmatrix}$. Finalmente, anularemos el elemento de la segunda fila y la segunda columna mendiante $\begin{matrix}f_1 \rightarrow f_1 \\f_2\rightarrow f_2\\-\frac{5}{4}f_2+f_3\rightarrow f_3\end{matrix}$ lo cual da lugar a la matriz reducida por Gauss $\begin{pmatrix}3&-1&1\\0&\frac{4}{3}&\frac{2}{3}\\0&0&-\frac{3}{2}\end{pmatrix}$, que, obviamente es una matriz triangular superior, y a la que denominaremos $U$.

Utilidades de la factorización $A=LU$

La factorización de una matriz regular como producto de una matriz triangular inferior y una matriz triangular superior facilita la resolución de un sistema de ecuaciones lineales, así como diversos cálculos matriciales. Veamos tres utilidades importantes:

• Dado un sistema de ecuaciones lineales $Ax=b$, tenemos que, al factorizar $A$ de este modo podemos escribir $LUx=b$, esto es $L(Ux)=b$, con lo cual podemos resolver el sistema en dos pasos: $\left\{\begin{matrix}\text{Primer paso: resolvemos}\quad Ly=b\\\text{Segundo paso: conocido}\quad y \quad \text{resolvemos:} \quad Ux=y \quad \text{y acabamos calculando} \quad x\end{matrix}\right.$
• Cálculo de la matriz inversa de $A$: Como $A=LU$, se tiene que $A^{-1}=(LU)^{-1}=U^{-1}\,L^{-1}$
• Cálculo del determinante de $A$: Como $A=LU$, y teniendo en cuenta que el determinante de un producto de matrices es igual al producto de sus determinantes se tiene que $\text{det}(A)=\text{det}(L)\cdot \text{det}(U)$. Entonces, si los elementos de la diagonal de $L$ son unos (factorización de Doolittle), resulta que $\text{det}(A)=1\cdot \displaystyle \prod_{i=1}^{n}\,u_{ii}=\displaystyle \prod_{i=1}^{n}\,u_{ii}$; y en el caso de que los unos estén en la diagonal de $U$ (factorización de Crout), tendremos que $\text{det}(A)=\displaystyle \prod_{i=1}^{n}\,\ell_{ii}\cdot 1=\displaystyle \prod_{i=1}^{n}\,\ell_{ii}$

Formalicemos un poco lo que acabamos de hacer acerca de cómo calcular la factorización $A=LU$:

Notemos ahora que, de los multiplicadores de la transformación, podemos escribir esta otra matriz triangular inferior y a la que denominaremos $L$ y escribiremos $L=\begin{pmatrix}1&0&0\\\frac{1}{3}&1&0\\\frac{2}{3}&\frac{5}{4}&1\end{pmatrix}$, donde hemos cambiado el signo de los números distintos de cero y de uno que hemos utilizado para efectuar las transformaciones. Pues bien, puede comprobarse que haciendo el producto $LU$ se obtiene la matriz $A$; en efecto $$\begin{pmatrix}1&0&0\\\frac{1}{3}&1&0\\\frac{2}{3}&\frac{5}{4}&1\end{pmatrix}\begin{pmatrix}3&-1&1\\0&\frac{4}{3}&\frac{2}{3}\\0&0&-\frac{3}{2}\end{pmatrix}=\begin{pmatrix}3&-1&1\\1&1&1\\2&1&0\end{pmatrix}$$ Esta es pues una manera sencilla de efectuar la factorización de una matriz regular dada $A$ como producto de una matriz triangular inferior por una matriz triangular superior, $A=LU$, siempre que no tengamos que utilizar técnicas de pivotamiento parcial —por encontrarnos con ceros en los elementos de la diagonal principal en las etapas del proceso de reducción, en cuyo caso deberemos extender esta técnica, tal y como ya trataré más adelante artícul—. Veremos también las utilidades de dicha factorización en el próximo artículo.

Cálculo de la matriz triangular inferior $L$

Démonos cuenta de que la matriz obtenida $E$ (la matriz de multiplicadores) puede a su vez descomponerse como un producto de dos matrices que corresponden a las dos etapas de que consta el proceso de reducción de esta matriz de orden $3$ y que vamos a denotar como $E(1)=\begin{pmatrix}1&0&0\\-\frac{1}{3}&1&0\\-\frac{2}{3}&0&1\end{pmatrix}$ y $E(2)=\begin{pmatrix}1&0&0\\0&1&0\\0& -\frac{5}{4}&1\end{pmatrix}$, y sus inversas, $E(1)^{-1}=\begin{pmatrix}1&0&0\\\frac{1}{3}&1&0\\\frac{2}{3}&0&1\end{pmatrix}$ y $E(2)^{-1}=\begin{pmatrix}1&0&0\\0&1&0\\0&\frac{5}{4}&1\end{pmatrix}$ ;donde, si recordemos lo que hemos hecho, y con vistas a formalizar un método para poder hacer descomposiciones de matrices de orden cualquiera, $E(1)=\begin{pmatrix}1&0&0\\-\frac{a_{21}^{(1)}}{a_{11}^{(1)}}&1&0\\-\frac{a_{31}^{(1)}}{a_{11}^{(1)}}&0&1\end{pmatrix}$ y $E(2)=\begin{pmatrix}1&0&0\\0&1&0\\0& -\frac{a_{32}^{(2)}}{a_{22}^{(2)}}&1\end{pmatrix}$, y sus inversas —se demuestra fácilmente que la matriz inversa $E(j)^{-1}$ se obtiene cambiando simplemente el signo de los elementos que corresponden a los multiplicadores de la matriz $E(j)$—, esto es, $E(1)^{-1}=\begin{pmatrix}1&0&0\\ \frac{a_{21}^{(1)}}{a_{11}^{(1)}}&1&0\\ \frac{a_{31}^{(1)}}{a_{11}^{(1)}}&0&1\end{pmatrix}$ y $E(2)^{-1}=\begin{pmatrix}1&0&0\\0&1&0\\0& \frac{a_{32}^{(2)}}{a_{22}^{(2)}}&1\end{pmatrix}$ (los superíndices —entre paréntesis— indican el orden de la etapa de reducción), con lo cual, $L=E(1)^{-1}\,E(2)^{-1}=\begin{pmatrix}1&0&0\\\frac{1}{3}&1&0\\\frac{2}{3}&0&1\end{pmatrix}\begin{pmatrix}1&0&0\\0&1&0\\0&\frac{5}{4}&1\end{pmatrix}=\begin{pmatrix}1&0&0\\\frac{1}{3}&1&0\\\frac{2}{3}&\frac{5}{4}&1\end{pmatrix}$.

Cálculo de la matriz triangular superior $U$

Observemos que, denotando por $A(1)\equiv A$, al ser $A$ una matriz de orden $3$, nos encontramos con que terminamos obteniendo la matriz triangular superior $U$ en dos etapas $\left\{\begin{matrix}A(2)=E(1)\,A(1)\\ U \equiv A(3)=E(2)\,A(2)\end{matrix}\right\}$, y terminamos. Si en lugar de simplemente $2$ etapas, tuviésemos $n-1$ (caso de una matriz de orden $n$), seguieremos la recurrencia $A(j+1)=E(j)\,A(j)$ para $j=1,2,\ldots,n-1$, siendo —recordémoslo— $A(1)\equiv A$ y $A(n)\equiv U$.

-oOo-

Observación: Conociendo $L$, para calcular $U$, fijémonos en que de la factorización que buscamos, $A=LU$, tenemos que, multiplicando ambos miembros por $L^{-1}$ por la izquierda, se tiene que $U=L^{-1}\,A=(E(1)^{-1}\,E(2)^{-1})^{-1}\,A=(E(2)^{-1})^{-1}\,(E(1)^{-1})^{-1}\,A=E(2)\,E(1)\,A$, con lo cual la matriz triangular superior $U$ es, como debe ser —valga como comprobación—, $$U=\begin{pmatrix}1&0&0\\0&1&0\\0&-\frac{5}{4}&1\end{pmatrix}\begin{pmatrix}1&0&0\\-\frac{1}{3}&1&0\\-\frac{2}{3}&0&1\end{pmatrix}\begin{pmatrix}3&-1&1\\1&1&1\\2&1&0\end{pmatrix}=\begin{pmatrix}3&-1&1\\0&\frac{4}{3}&\frac{2}{3}\\0&0&-\frac{3}{2}\end{pmatrix}$$

-oOo-

Así, en el caso de factorizar una matriz de orden $4$, la matriz triangular inferior haremos el siguiente cálculo, $L=E(1)\,E(2)\,E(3)$; y para calcular la matriz triangular superior $U$ tendríamos tres pasos en lugar de dos: $\left\{\begin{matrix} A(2)=E(1)\,A(1)\\ A(3)=E(2)\,A(2) \\ U\equiv A(4)=E(3)\,A(3) \end{matrix}\right\}$ donde $A(1) \equiv A$, con las respectivas matrices de multiplicadores (ahora hay tres): $E(1)=\begin{pmatrix}1&0&0&0 \\ -\frac{a_{21}^{(1)}}{a_{11}^{(1)}} & 1 & 0 & 0 \\ -\frac{a_{31}^{(1)}}{a_{11}^{(1)}} & 0 & 1 & 0 \\ -\frac{a_{41}^{(1)}}{a_{11}^{(1)}} & 0 & 0 & 1 \end{pmatrix}$, $E(2)=\begin{pmatrix}1& 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & -\frac{a_{32}^{(2)}}{a_{22}^{(2)}} & 1 & 0 \\ 0 & -\frac{a_{42}^{(2)}}{a_{22}^{(2)}} & 0 & 1 \end{pmatrix}$ y $E(3)=\begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & -\frac{a_{43}^{(3)}}{a_{33}^{(3)}} & 1\end{pmatrix}$, así como sus respectivas inversas: $E(1)^{-1}=\begin{pmatrix}1&0&0&0 \\ \frac{a_{21}^{(1)}}{a_{11}^{(1)}} & 1 & 0 & 0 \\ \frac{a_{31}^{(1)}}{a_{11}^{(1)}} & 0 & 1 & 0 \\ \frac{a_{41}^{(1)}}{a_{11}^{(1)}} & 0 & 0 & 1 \end{pmatrix}$, $E(2)^{-1}=\begin{pmatrix}1& 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & \frac{a_{32}^{(2)}}{a_{22}^{(2)}} & 1 & 0 \\ 0 & \frac{a_{42}^{(2)}}{a_{22}^{(2)}} & 0 & 1 \end{pmatrix}$ y $E(3)^{-1}=\begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & \frac{a_{43}^{(3)}}{a_{33}^{(3)}} & 1\end{pmatrix}$

Generalización. Factorización gaussiana $LU$ de una matriz de orden $n$ (sin que sea necesario el pivotaje)

A partir de lo que se acaba de decir, resulta ya muy fáoil generalizar lo que se ha encontrado para una matriz $A$ de orden $n$ —siempre que no haya que recurrir a permutar filas (pivotamiento) por encontrarnos con ceros o bien con elementos próximos a cero en las posiciones diagonales—, podemos escribir:

1. $L=E(1)^{-1}\,E(2)^{-1}\,\ldots\,E(n-1)^{-1}$
2. $U\equiv A(n)=E(n-1)\,A(n-1)=E(n-1)\,E(n-2)\,A(n-2)$
   $=\ldots=E(n-1)\,E(n-2)\,\ldots\,E(2)\,\,E(1)\,A(1)$, donde $A(1)\equiv A$

donde $$E(j)=\begin{pmatrix} 1&0&\ldots &0&\ldots &0&0&0\\ 0&1&\ldots &0&\ldots &0&0&0\\ \vdots & \vdots & \ddots &\vdots &\vdots &\vdots &\vdots &\vdots\\ \vdots & \vdots & \vdots &1 &\vdots &\vdots &\vdots &\vdots\\ 0 & 0 & \ldots &-\frac{a_{j+1,j}^{(j)}}{a_{jj}^{(j)}} &1 &\ldots &\ldots &0\\ \vdots & \vdots & \vdots &\vdots &\vdots &\ddots &\vdots &\vdots\\ 0 & 0 & \ldots &-\frac{a_{n-1,j}^{(j)}}{a_{jj}^{(j)}} &\ldots &\ldots &1 &0\\ 0 & 0 & \ldots &-\frac{a_{n,j}^{(j)}}{a_{jj}^{(j)}} &\ldots &\ldots &0 &1\\ \end{pmatrix}$$ y su respectiva matriz inversa $$E(j)^{-1}=\begin{pmatrix} 1&0&\ldots &0&\ldots &0&0&0\\ 0&1&\ldots &0&\ldots &0&0&0\\ \vdots & \vdots & \ddots &\vdots &\vdots &\vdots &\vdots &\vdots\\ \vdots & \vdots & \vdots &1 &\vdots &\vdots &\vdots &\vdots\\ 0 & 0 & \ldots &\frac{a_{j+1,j}^{(j)}}{a_{jj}^{(j)}} &1 &\ldots &\ldots &0\\ \vdots & \vdots & \vdots &\vdots &\vdots &\ddots &\vdots &\vdots\\ 0 & 0 & \ldots &\frac{a_{n-1,j}^{(j)}}{a_{jj}^{(j)}} &\ldots &\ldots &1 &0\\ 0 & 0 & \ldots &\frac{a_{n,j}^{(j)}}{a_{jj}^{(j)}} &\ldots &\ldots &0 &1\\ \end{pmatrix}$$ de donde se obtiene la matriz triangular inferior $$L=\begin{pmatrix} 1&0&\ldots & 0 &\ldots & 0 & 0 & 0 \\ \frac{a_{21}^{(1)}}{a_{11}^{(1)}}&1&\ldots & 0 &\ldots & 0 & 0 & 0 \\ \vdots&\vdots&\ddots & \vdots &\vdots & \vdots & \vdots & \vdots \\ \vdots&\vdots&\vdots & 1 &\vdots & \vdots & \vdots & \vdots \\ \frac{a_{j1}^{(1)}}{a_{11}^{(1)}}& \frac{a_{j2}^{(2)}}{a_{22}^{(2)}}&\ldots & \frac{a_{j+1,j}^{(j)}}{a_{jj}^{(j)}} & 1 & \ldots & \ldots & 0 \\ \vdots&\vdots&\vdots & \vdots &\vdots & \ddots & \vdots & \vdots \\ \frac{a_{n-1,1}^{(1)}}{a_{11}^{(1)}}& \frac{a_{n-1,2}^{(2)}}{a_{22}^{(2)}}&\ldots & \frac{a_{n-1,j}^{(j)}}{a_{jj}^{(j)}} & \ldots & \ldots & 1 & 0 \\ \frac{a_{n,1}^{(1)}}{a_{11}^{(1)}}& \frac{a_{n,2}^{(2)}}{a_{22}^{(2)}}&\ldots & \frac{a_{n,j}^{(j)}}{a_{jj}^{(j)}} & \ldots & \ldots & \frac{a_{n,n-1}^{(n-1)}}{a_{n-1,n-1}^{(n-1)}} & 1 \\ \end{pmatrix}$$ $\diamond$

-oOo-

Referencias:
[1] A. Aubanell, A. Benseny y A. Deslshams, Útiles de cálculo numérico (Labor, Barcelona, 1993).
[2] C. Moreno, Introducción al cálculo numérico (UNED, Madrid, 2011).