Descomposición LU

Tipo de factorización matricial

En el análisis numérico y el álgebra lineal , la descomposición LU o factorización factoriza una matriz como el producto de una matriz triangular inferior y una matriz triangular superior (ver descomposición matricial ). El producto a veces también incluye una matriz de permutación . La descomposición LU puede verse como la forma matricial de la eliminación gaussiana . Las computadoras generalmente resuelven sistemas cuadrados de ecuaciones lineales utilizando la descomposición LU, y también es un paso clave al invertir una matriz o calcular el determinante de una matriz. La descomposición LU fue introducida por el astrónomo polaco Tadeusz Banachiewicz en 1938. ^[1] Para citar: "Parece que Gauss y Doolittle aplicaron el método [de eliminación] solo a ecuaciones simétricas. Autores más recientes, por ejemplo, Aitken, Banachiewicz, Dwyer y Crout ... han enfatizado el uso del método, o variaciones del mismo, en conexión con problemas no simétricos ... Banachiewicz ... vio el punto ... de que el problema básico es realmente uno de factorización de matrices, o "descomposición", como él lo llamó". ^[2] También se lo conoce a veces como descomposición LR (factores en matrices triangulares izquierdas y derechas).

Definiciones

Descomposición LDU de una matriz de Walsh

Sea A una matriz cuadrada. Una factorización LU se refiere a la factorización de A , con ordenamientos o permutaciones de filas y/o columnas adecuados, en dos factores: una matriz triangular inferior L y una matriz triangular superior U :

${\displaystyle A=LU.}$

En la matriz triangular inferior, todos los elementos que se encuentran por encima de la diagonal son cero; en la matriz triangular superior, todos los elementos que se encuentran por debajo de la diagonal son cero. Por ejemplo, para una matriz A de 3 × 3 , su descomposición LU se ve así:

${\displaystyle {\begin{bmatrix}a_{11}&a_{12}&a_{13}\\a_{21}&a_{22}&a_{23}\\a_{31}&a_{32}&a_{33}\end{bmatrix}}={\begin{bmatrix}\ell _{11}&0&0\\\ell _{21}&\ell _{22}&0\\\ell _{31}&\ell _{32}&\ell _{33}\end{bmatrix}}{\begin{bmatrix}u_{11}&u_{12}&u_{13}\\0&u_{22}&u_{23}\\0&0&u_{33}\end{bmatrix}}.}$

Sin un ordenamiento adecuado o permutaciones en la matriz, la factorización puede no materializarse. Por ejemplo, es fácil verificar (al expandir la multiplicación de matrices ) que . Si , entonces al menos uno de y tiene que ser cero, lo que implica que L o U es singular . Esto es imposible si A no es singular (invertible). Este es un problema de procedimiento. Puede eliminarse simplemente reordenando las filas de A de modo que el primer elemento de la matriz permutada no sea cero. El mismo problema en los pasos de factorización posteriores puede eliminarse de la misma manera; consulte el procedimiento básico a continuación. ${\textstyle a_{11}=\ell _{11}u_{11}}$ ${\textstyle a_{11}=0}$ ${\textstyle \ell _{11}}$ ${\textstyle u_{11}}$

Factorización LU con pivoteo parcial

Resulta que una permutación adecuada en filas (o columnas) es suficiente para la factorización LU. La factorización LU con pivoteo parcial (LUP) se refiere a menudo a la factorización LU con permutaciones de filas únicamente:

${\displaystyle PA=LU,}$

donde L y U son nuevamente matrices triangulares inferior y superior, y P es una matriz de permutación que, cuando se multiplica por la izquierda a A , reordena las filas de A. Resulta que todas las matrices cuadradas se pueden factorizar de esta forma, ^[3] y la factorización es numéricamente estable en la práctica. ^[4] Esto hace que la descomposición LUP sea una técnica útil en la práctica.

Factorización LU con pivoteo completo

Una factorización LU con pivoteo completo implica permutaciones tanto de filas como de columnas:

${\displaystyle PAQ=LU,}$

donde L , U y P se definen como antes, y Q es una matriz de permutación que reordena las columnas de A. ^[5 ]

Descomposición diagonal inferior superior (LDU)

Una descomposición diagonal inferior superior (LDU) es una descomposición de la forma

${\displaystyle A=LDU,}$

donde D es una matriz diagonal , y L y U son matrices triangulares unitarias , lo que significa que todas las entradas en las diagonales de L y U son uno.

Matrices rectangulares

Anteriormente, exigimos que A sea una matriz cuadrada, pero todas estas descomposiciones se pueden generalizar también a matrices rectangulares. ^[6] En ese caso, L y D son matrices cuadradas, ambas con el mismo número de filas que A , y U tiene exactamente las mismas dimensiones que A . La forma triangular superior debe interpretarse como que tiene solo cero entradas debajo de la diagonal principal, que comienza en la esquina superior izquierda. De manera similar, el término más preciso para U es que es la forma escalonada por filas de la matriz A .

Ejemplo

Factorizamos la siguiente matriz de 2 por 2:

${\displaystyle {\begin{bmatrix}4&3\\6&3\end{bmatrix}}={\begin{bmatrix}\ell _{11}&0\\\ell _{21}&\ell _{22}\end{bmatrix}}{\begin{bmatrix}u_{11}&u_{12}\\0&u_{22}\end{bmatrix}}.}$

Una forma de encontrar la descomposición LU de esta matriz simple sería simplemente resolver las ecuaciones lineales por inspección. Al expandir la multiplicación de matrices se obtiene

${\displaystyle {\begin{aligned}\ell _{11}\cdot u_{11}+0\cdot 0&=4\\\ell _{11}\cdot u_{12}+0\cdot u_{22}&=3\\\ell _{21}\cdot u_{11}+\ell _{22}\cdot 0&=6\\\ell _{21}\cdot u_{12}+\ell _{22}\cdot u_{22}&=3.\end{aligned}}}$

Este sistema de ecuaciones está indeterminado . En este caso, dos elementos cualesquiera no nulos de las matrices L y U son parámetros de la solución y pueden asignarse arbitrariamente a cualquier valor no nulo. Por lo tanto, para encontrar la descomposición LU única, es necesario poner alguna restricción a las matrices L y U. Por ejemplo, podemos exigir convenientemente que la matriz triangular inferior L sea una matriz triangular unitaria, de modo que todas las entradas de su diagonal principal se establezcan en uno. Entonces, el sistema de ecuaciones tiene la siguiente solución:

${\displaystyle {\begin{aligned}\ell _{11}=\ell _{22}&=1\\\ell _{21}&=1.5\\u_{11}&=4\\u_{12}&=3\\u_{22}&=-1.5\end{aligned}}}$

Sustituyendo estos valores en la descomposición LU anterior se obtiene

${\displaystyle {\begin{bmatrix}4&3\\6&3\end{bmatrix}}={\begin{bmatrix}1&0\\1.5&1\end{bmatrix}}{\begin{bmatrix}4&3\\0&-1.5\end{bmatrix}}.}$

Existencia y singularidad

Matrices cuadradas

Toda matriz cuadrada admite factorizaciones LUP y PLU . ^[3] Si es invertible , entonces admite una factorización LU (o LDU ) si y solo si todos sus principales menores ^[7] son ​​distintos de cero ^[8] (por ejemplo no admite una factorización LU o LDU ). Si es una matriz singular de rango , entonces admite una factorización LU si los primeros principales menores son distintos de cero, aunque lo inverso no es cierto. ^[9] ${\textstyle A}$ ${\textstyle A}$ ${\displaystyle {\begin{bmatrix}0&1\\1&0\end{bmatrix}}}$ ${\textstyle A}$ ${\textstyle k}$ ${\textstyle k}$

Si una matriz cuadrada invertible tiene una LDU (factorización con todas las entradas diagonales de L y U iguales a 1), entonces la factorización es única. ^[8] En ese caso, la factorización LU también es única si requerimos que la diagonal de (o ) consista en unos. ${\textstyle L}$ ${\textstyle U}$

En general, cualquier matriz cuadrada podría tener uno de los siguientes: ${\displaystyle A_{n\times n}}$

1. una factorización LU única (como se mencionó anteriormente);
2. infinitas factorizaciones LU si dos o más de las primeras ( n −1) columnas son linealmente dependientes o cualquiera de las primeras ( n −1) columnas son 0;
3. no hay factorización LU si las primeras ( n −1) columnas no son cero y son linealmente independientes y al menos un menor principal es cero.

En el caso 3, se puede aproximar una factorización LU cambiando una entrada diagonal para evitar un menor principal inicial cero. ^[10] ${\displaystyle a_{jj}}$ ${\displaystyle a_{jj}\pm \varepsilon }$

Matrices simétricas definidas positivas

Si A es una matriz simétrica (o hermítica , si A es compleja) definida positivamente , podemos ordenar las cosas de modo que U sea la transpuesta conjugada de L. Es decir, podemos escribir A como

${\displaystyle A=LL^{*}.\,}$

Esta descomposición se denomina descomposición de Cholesky . Si es definida positiva, entonces la descomposición de Cholesky existe y es única. Además, el cálculo de la descomposición de Cholesky es más eficiente y numéricamente más estable que el cálculo de otras descomposiciones LU. ${\displaystyle A}$

Matrices generales

Para una matriz (no necesariamente invertible) sobre cualquier cuerpo, se conocen las condiciones exactas necesarias y suficientes bajo las cuales tiene una factorización LU. Las condiciones se expresan en términos de los rangos de ciertas submatrices. El algoritmo de eliminación gaussiana para obtener la descomposición LU también se ha extendido a este caso más general. ^[11]

Algoritmos

Fórmula cerrada

Cuando existe una factorización LDU y es única, hay una fórmula cerrada (explícita) para los elementos de L , D y U en términos de razones de determinantes de ciertas submatrices de la matriz original A . ^[12] En particular, , y para , es la razón de la -ésima submatriz principal a la -ésima submatriz principal. El cálculo de los determinantes es computacionalmente costoso , por lo que esta fórmula explícita no se utiliza en la práctica. ${\textstyle D_{1}=A_{1,1}}$ ${\textstyle i=2,\ldots ,n}$ ${\textstyle D_{i}}$ ${\textstyle i}$ ${\textstyle (i-1)}$

Utilizando la eliminación gaussiana

El siguiente algoritmo es esencialmente una forma modificada de eliminación gaussiana . Calcular una descomposición LU utilizando este algoritmo requiere operaciones de punto flotante, ignorando los términos de orden inferior. El pivoteo parcial agrega solo un término cuadrático; este no es el caso del pivoteo completo. ^[13] ${\displaystyle {\tfrac {2}{3}}n^{3}}$

Explicación generalizada

Notación

Dada una matriz N × N , se define como la versión original, sin modificar, de la matriz . El superíndice entre paréntesis (por ejemplo, ) de la matriz es la versión de la matriz. La matriz es la matriz en la que los elementos debajo de la diagonal principal ya se han eliminado a 0 mediante eliminación gaussiana para las primeras columnas. ${\displaystyle A=(a_{i,j})_{1\leq i,j\leq N}}$ ${\displaystyle A^{(0)}}$ ${\displaystyle A}$ ${\displaystyle (0)}$ ${\displaystyle A}$ ${\displaystyle A^{(n)}}$ ${\displaystyle A}$ ${\displaystyle n}$

A continuación se muestra una matriz para observar que nos ayudará a recordar la notación (donde cada uno representa cualquier número real en la matriz): ${\displaystyle *}$

${\displaystyle A^{(n-1)}={\begin{pmatrix}*&&&\cdots &&&*\\0&\ddots &&&&\\&\ddots &*&&&\\\vdots &&0&a_{n,n}^{(n-1)}&&&\vdots \\&&\vdots &a_{i,n}^{(n-1)}&*\\&&&\vdots &\vdots &\ddots \\0&\cdots &0&a_{i,n}^{(n-1)}&*&\cdots &*\end{pmatrix}}}$

Procedimiento

Durante este proceso, modificamos gradualmente la matriz mediante operaciones de fila hasta que se convierte en la matriz en la que todos los elementos debajo de la diagonal principal son iguales a cero. Durante este proceso, crearemos simultáneamente dos matrices separadas y , de modo que . ${\displaystyle A}$ ${\displaystyle U}$ ${\displaystyle P}$ ${\displaystyle L}$ ${\displaystyle PA=LU}$

Definimos la matriz de permutación final como la matriz identidad que tiene todas las mismas filas intercambiadas en el mismo orden que la matriz mientras se transforma en la matriz . Para nuestra matriz , podemos comenzar intercambiando filas para proporcionar las condiciones deseadas para la n-ésima columna. Por ejemplo, podemos intercambiar filas para realizar un pivoteo parcial, o podemos hacerlo para establecer el elemento pivote en la diagonal principal en un número distinto de cero para que podamos completar la eliminación gaussiana. ${\displaystyle P}$ ${\displaystyle A}$ ${\displaystyle U}$ ${\displaystyle A^{(n-1)}}$ ${\displaystyle a_{n,n}}$

Para nuestra matriz , queremos establecer cada elemento a continuación en cero (donde es el elemento en la n-ésima columna de la diagonal principal). Denotaremos cada elemento a continuación como (donde ). Para establecer en cero, establecemos para cada fila . Para esta operación, . Una vez que hemos realizado las operaciones de fila para las primeras columnas, hemos obtenido una matriz triangular superior que se denota por . ${\displaystyle A^{(n-1)}}$ ${\displaystyle a_{n,n}^{(n-1)}}$ ${\displaystyle a_{n,n}^{(n-1)}}$ ${\displaystyle a_{n,n}^{(n-1)}}$ ${\displaystyle a_{i,n}^{(n-1)}}$ ${\displaystyle i=n+1,\dotsc ,N}$ ${\displaystyle a_{i,n}^{(n-1)}}$ ${\displaystyle row_{i}=row_{i}-(\ell _{i,n})\cdot row_{n}}$ ${\displaystyle i}$ ${\textstyle \ell _{i,n}:={a_{i,n}^{(n-1)}}/{a_{n,n}^{(n-1)}}}$ ${\displaystyle N-1}$ ${\displaystyle A^{(N-1)}}$ ${\displaystyle U}$

También podemos crear la matriz triangular inferior denotada como , ingresando directamente los valores previamente calculados de mediante la fórmula siguiente. ${\textstyle L}$ ${\displaystyle \ell _{i,n}}$

${\displaystyle L={\begin{pmatrix}1&0&\cdots &0\\\ell _{2,1}&\ddots &\ddots &\vdots \\\vdots &\ddots &\ddots &0\\\ell _{N,1}&\cdots &\ell _{N,N-1}&1\end{pmatrix}}}$

Ejemplo

Si nos dan la matriz, elegiremos implementar pivoteo parcial y, por lo tanto, intercambiar la primera y la segunda fila para que nuestra matriz y la primera iteración de nuestra matriz se conviertan respectivamente en Una vez que hemos intercambiado las filas, podemos eliminar los elementos debajo de la diagonal principal en la primera columna realizando de modo que, Una vez que se han restado estas filas, hemos derivado de la matriz Debido a que estamos implementando pivoteo parcial, intercambiamos la segunda y la tercera filas de nuestra matriz derivada y la versión actual de nuestra matriz respectivamente para obtener Ahora, eliminamos los elementos debajo de la diagonal principal en la segunda columna realizando de modo que . Debido a que no existen elementos distintos de cero debajo de la diagonal principal en nuestra iteración actual de después de esta resta de filas, esta resta de filas deriva nuestra matriz final (denotada como ) y matriz final: Después de intercambiar también las filas correspondientes, obtenemos nuestra matriz final: Ahora estas matrices tienen una relación tal que . $${\displaystyle A={\begin{pmatrix}0&5&{\frac {22}{3}}\\4&2&1\\2&7&9\\\end{pmatrix}},}$$ ${\displaystyle A}$ ${\displaystyle P}$ $${\displaystyle A^{(0)}={\begin{pmatrix}4&2&1\\0&5&{\frac {22}{3}}\\2&7&9\\\end{pmatrix}},\quad P^{(0)}={\begin{pmatrix}0&1&0\\1&0&0\\0&0&1\\\end{pmatrix}}.}$$ $${\displaystyle {\begin{alignedat}{0}row_{2}=row_{2}-(\ell _{2,1})\cdot row_{1}\\row_{3}=row_{3}-(\ell _{3,1})\cdot row_{1}\end{alignedat}}}$$ $${\displaystyle {\begin{alignedat}{0}\ell _{2,1}={\frac {0}{4}}=0\\\ell _{3,1}={\frac {2}{4}}=0.5\end{alignedat}}}$$ ${\displaystyle A^{(1)}}$ $${\displaystyle A^{(1)}={\begin{pmatrix}4&2&1\\0&5&{\frac {22}{3}}\\0&6&8.5\\\end{pmatrix}}.}$$ ${\displaystyle P}$ $${\displaystyle A^{(1)}={\begin{pmatrix}4&2&1\\0&6&8.5\\0&5&{\frac {22}{3}}\\\end{pmatrix}},\quad P^{(1)}={\begin{pmatrix}0&1&0\\0&0&1\\1&0&0\\\end{pmatrix}}.}$$ ${\displaystyle row_{3}=row_{3}-(\ell _{3,2})\cdot row_{2}}$ ${\textstyle \ell _{3,2}={\frac {5}{6}}}$ ${\displaystyle A}$ ${\displaystyle A}$ ${\displaystyle U}$ ${\displaystyle P}$ $${\displaystyle A^{(2)}=A^{(N-1)}=U={\begin{pmatrix}4&2&1\\0&6&8.5\\0&0&0.25\\\end{pmatrix}},\quad P={\begin{pmatrix}0&1&0\\0&0&1\\1&0&0\\\end{pmatrix}}.}$$ ${\displaystyle L}$ $${\displaystyle L={\begin{pmatrix}1&0&0\\\ell _{3,1}&1&0\\\ell _{2,1}&\ell _{3,2}&1\\\end{pmatrix}}={\begin{pmatrix}1&0&0\\0.5&1&0\\0&{\frac {5}{6}}&1\\\end{pmatrix}}}$$ ${\displaystyle PA=LU}$

Relaciones cuando no se intercambian filas

Si no intercambiamos filas en absoluto durante este proceso, podemos realizar las operaciones de fila simultáneamente para cada columna estableciendo donde es la matriz identidad N × N con su n -ésima columna reemplazada por el vector transpuesto En otras palabras, la matriz triangular inferior ${\displaystyle n}$ ${\displaystyle A^{(n)}:=L_{n}^{-1}A^{(n-1)},}$ ${\displaystyle L_{n}^{-1}}$ ${\displaystyle {\begin{pmatrix}0&\dotsm &0&1&-\ell _{n+1,n}&\dotsm &-\ell _{N,n}\end{pmatrix}}^{\textsf {T}}.}$

${\displaystyle L_{n}^{-1}={\begin{pmatrix}1&&&&&\\&\ddots &&&&\\&&1&&&\\&&-\ell _{n+1,n}&&&\\&&\vdots &&\ddots &\\&&-\ell _{N,n}&&&1\end{pmatrix}}.}$

Realizar todas las operaciones de fila para las primeras columnas utilizando la fórmula es equivalente a encontrar la descomposición Denote de modo que . ${\displaystyle N-1}$ ${\displaystyle A^{(n)}:=L_{n}^{-1}A^{(n-1)}}$ $${\displaystyle A=L_{1}L_{1}^{-1}A^{(0)}=L_{1}A^{(1)}=L_{1}L_{2}L_{2}^{-1}A^{(1)}=L_{1}L_{2}A^{(2)}=\dotsm =L_{1}\dotsm L_{N-1}A^{(N-1)}.}$$ ${\textstyle L=L_{1}\dotsm L_{N-1}}$ ${\displaystyle A=LA^{(N-1)}=LU}$

Ahora calculemos la secuencia de . Sabemos que tiene la siguiente fórmula. ${\displaystyle L_{1}\dotsm L_{N-1}}$ ${\displaystyle L_{i}}$

${\displaystyle L_{n}={\begin{pmatrix}1&&&&&\\&\ddots &&&&\\&&1&&&\\&&\ell _{n+1,n}&&&\\&&\vdots &&\ddots &\\&&\ell _{N,n}&&&1\end{pmatrix}}}$

Si hay dos matrices triangulares inferiores con 1 en la diagonal principal y ninguna tiene un elemento distinto de cero debajo de la diagonal principal en la misma columna que la otra, entonces podemos incluir todos los elementos distintos de cero en su misma ubicación en el producto de las dos matrices. Por ejemplo:

${\displaystyle \left({\begin{array}{ccccc}1&0&0&0&0\\77&1&0&0&0\\12&0&1&0&0\\63&0&0&1&0\\7&0&0&0&1\end{array}}\right)\left({\begin{array}{ccccc}1&0&0&0&0\\0&1&0&0&0\\0&22&1&0&0\\0&33&0&1&0\\0&44&0&0&1\end{array}}\right)=\left({\begin{array}{ccccc}1&0&0&0&0\\77&1&0&0&0\\12&22&1&0&0\\63&33&0&1&0\\7&44&0&0&1\end{array}}\right)}$

Finalmente, multiplicamos y generamos la matriz fusionada denotada como (como se mencionó anteriormente). Usando la matriz , obtenemos ${\displaystyle L_{i}}$ ${\textstyle L}$ ${\textstyle L}$ ${\displaystyle A=LU.}$

Está claro que para que este algoritmo funcione, es necesario tener en cada paso (ver la definición de ). Si esta suposición falla en algún punto, es necesario intercambiar la fila n -ésima con otra fila debajo de ella antes de continuar. Es por esto que una descomposición LU en general se ve como . ${\displaystyle a_{n,n}^{(n-1)}\neq 0}$ ${\displaystyle \ell _{i,n}}$ ${\displaystyle P^{-1}A=LU}$

Descomposición de LU Crout

Nótese que la descomposición obtenida mediante este procedimiento es una descomposición de Doolittle : la diagonal principal de L está compuesta únicamente por 1 s. Si se procede eliminando elementos por encima de la diagonal principal agregando múltiplos de las columnas (en lugar de eliminar elementos por debajo de la diagonal agregando múltiplos de las filas ), obtendremos una descomposición de Crout , donde la diagonal principal de U es de 1 s.

Otra forma (equivalente) de producir una descomposición de Crout de una matriz dada A es obtener una descomposición de Doolittle de la transpuesta de A . De hecho, si es la descomposición LU obtenida a través del algoritmo presentado en esta sección, entonces al tomar y , tenemos que es una descomposición de Crout. ${\textstyle A^{\textsf {T}}=L_{0}U_{0}}$ ${\textstyle L=U_{0}^{\textsf {T}}}$ ${\textstyle U=L_{0}^{\textsf {T}}}$ ${\displaystyle A=LU}$

A través de la recursión

Cormen et al. ^[14] describen un algoritmo recursivo para la descomposición LUP.

Dada una matriz A , sea P ₁ una matriz de permutación tal que

${\displaystyle P_{1}A=\left({\begin{array}{c|ccc}a&&w^{\textsf {T}}&\\\hline &&&\\v&&A'&\\&&&\end{array}}\right)}$,

donde , si hay una entrada distinta de cero en la primera columna de A ; o tome P ₁ como la matriz identidad en caso contrario. Ahora sea , si ; o en caso contrario. Tenemos ${\textstyle a\neq 0}$ ${\textstyle c=1/a}$ ${\textstyle a\neq 0}$ ${\textstyle c=0}$

${\displaystyle P_{1}A=\left({\begin{array}{c|ccc}1&&0&\\\hline &&&\\cv&&I_{n-1}&\\&&&\end{array}}\right)\left({\begin{array}{c|c}a&w^{\textsf {T}}\\\hline &\\0&A'-cvw^{\textsf {T}}\\&\end{array}}\right).}$

Ahora podemos encontrar recursivamente una descomposición LUP . Sea . Por lo tanto ${\textstyle P'\left(A'-cvw^{\textsf {T}}\right)=L'U'}$ ${\textstyle v'=P'v}$

${\displaystyle \left({\begin{array}{c|ccc}1&&0&\\\hline &&&\\0&&P'&\\&&&\end{array}}\right)P_{1}A=\left({\begin{array}{c|ccc}1&&0&\\\hline &&&\\cv'&&L'&\\&&&\end{array}}\right)\left({\begin{array}{c|ccc}a&&w^{\textsf {T}}&\\\hline &&&\\0&&U'&\\&&&\end{array}}\right),}$

que es una descomposición LUP de A.

Algoritmo aleatorio

Es posible encontrar una aproximación de rango bajo a una descomposición LU utilizando un algoritmo aleatorio. Dada una matriz de entrada y un rango bajo deseado , la LU aleatoria devuelve matrices de permutación y matrices trapezoidales inferior/superior de tamaño y respectivamente, tales que con alta probabilidad , donde es una constante que depende de los parámetros del algoritmo y es el -ésimo valor singular de la matriz de entrada . ^[15] ${\textstyle A}$ ${\textstyle k}$ ${\textstyle P,Q}$ ${\textstyle L,U}$ ${\textstyle m\times k}$ ${\textstyle k\times n}$ ${\textstyle \left\|PAQ-LU\right\|_{2}\leq C\sigma _{k+1}}$ ${\textstyle C}$ ${\textstyle \sigma _{k+1}}$ ${\textstyle (k+1)}$ ${\textstyle A}$

Complejidad teórica

Si dos matrices de orden n se pueden multiplicar en el tiempo M ( n ), donde M ( n ) ≥ n ^a para algún a > 2, entonces se puede calcular una descomposición LU en el tiempo O( M ( n )). ^[16] Esto significa, por ejemplo, que existe un algoritmo O( n ^2.376 ) basado en el algoritmo Coppersmith–Winograd .

Descomposición en matriz dispersa

Se han desarrollado algoritmos especiales para factorizar matrices dispersas de gran tamaño . Estos algoritmos intentan encontrar los factores dispersos L y U. Idealmente, el costo del cálculo está determinado por la cantidad de entradas distintas de cero, en lugar de por el tamaño de la matriz.

Estos algoritmos utilizan la libertad de intercambiar filas y columnas para minimizar el relleno (entradas que cambian de un valor inicial cero a un valor distinto de cero durante la ejecución de un algoritmo).

El tratamiento general de los ordenamientos que minimizan el relleno se puede abordar utilizando la teoría de grafos .

Aplicaciones

Resolver ecuaciones lineales

Dado un sistema de ecuaciones lineales en forma matricial

${\displaystyle A\mathbf {x} =\mathbf {b} ,}$

Queremos resolver la ecuación para x , dados A y b . Supongamos que ya hemos obtenido la descomposición LUP de A tal que , entonces . ${\textstyle PA=LU}$ ${\textstyle LU\mathbf {x} =P\mathbf {b} }$

En este caso la solución se realiza en dos pasos lógicos:

1. Primero, resolvemos la ecuación para y . ${\textstyle L\mathbf {y} =P\mathbf {b} }$
2. En segundo lugar, resolvemos la ecuación para x . ${\textstyle U\mathbf {x} =\mathbf {y} }$

En ambos casos se trata de matrices triangulares ( L y U ), que pueden resolverse directamente mediante sustitución hacia adelante y hacia atrás sin utilizar el proceso de eliminación gaussiana (sin embargo, necesitamos este proceso o equivalente para calcular la descomposición LU en sí).

El procedimiento anterior se puede aplicar repetidamente para resolver la ecuación varias veces para diferentes b . En este caso, es más rápido (y más conveniente) hacer una descomposición LU de la matriz A una vez y luego resolver las matrices triangulares para los diferentes b , en lugar de usar la eliminación gaussiana cada vez. Se podría pensar que las matrices L y U han "codificado" el proceso de eliminación gaussiana.

El costo de resolver un sistema de ecuaciones lineales es aproximadamente el de operaciones de punto flotante si la matriz tiene un tamaño de . Esto hace que sea el doble de rápido que los algoritmos basados ​​en la descomposición QR , que cuesta aproximadamente el de operaciones de punto flotante cuando se utilizan reflexiones de Householder . Por esta razón, la descomposición LU suele ser la preferida. ^[17] ${\textstyle {\frac {2}{3}}n^{3}}$ ${\textstyle A}$ ${\textstyle n}$ ${\textstyle {\frac {4}{3}}n^{3}}$

Invertir una matriz

Al resolver sistemas de ecuaciones, b se suele considerar un vector con una longitud igual a la altura de la matriz A. Sin embargo, en la inversión de matrices, en lugar del vector b , tenemos la matriz B , donde B es una matriz n por p , de modo que estamos tratando de encontrar una matriz X (también una matriz n por p ):

${\displaystyle AX=LUX=B.}$

Podemos utilizar el mismo algoritmo presentado anteriormente para resolver cada columna de la matriz X. Ahora supongamos que B es la matriz identidad de tamaño n . Se deduciría que el resultado X debe ser la inversa de A. ^[18]

Calculando el determinante

Dada la descomposición LUP de una matriz cuadrada A , el determinante de A se puede calcular directamente como ${\textstyle A=P^{-1}LU}$

${\displaystyle \det(A)=\det \left(P^{-1}\right)\det(L)\det(U)=(-1)^{S}\left(\prod _{i=1}^{n}l_{ii}\right)\left(\prod _{i=1}^{n}u_{ii}\right).}$

La segunda ecuación se deriva del hecho de que el determinante de una matriz triangular es simplemente el producto de sus entradas diagonales, y que el determinante de una matriz de permutación es igual a (−1) ^S, donde S es el número de intercambios de filas en la descomposición.

En el caso de la descomposición LU con pivoteo completo, también es igual al lado derecho de la ecuación anterior, si dejamos que S sea el número total de intercambios de filas y columnas. ${\textstyle \det(A)}$

El mismo método se aplica fácilmente a la descomposición LU estableciendo P igual a la matriz identidad.

Ejemplos de código

Ejemplo de código C

/* ENTRADA: A - matriz de punteros a filas de una matriz cuadrada que tiene dimensión N * Tol - pequeño número de tolerancia para detectar fallas cuando la matriz está cerca de degenerarse * SALIDA: La matriz A se modifica, contiene una copia de ambas matrices LE y U como A=(LE)+U tal que P*A=L*U. * La matriz de permutación no se almacena como una matriz, sino en un vector entero P de tamaño N+1 * que contiene índices de columna donde la matriz de permutación tiene "1". El último elemento P[N]=S+N, * donde S es el número de intercambios de filas necesarios para el cálculo del determinante, det(P)=(-1)^S */ int LUPDecompose ( double ** A , int N , double Tol , int * P ) {          int i , j , k , imax ; doble maxA , * ptr , absA ;          para ( i = 0 ; i <= N ; i ++ ) P [ i ] = i ; // Matriz de permutación unitaria, P[N] inicializada con N            para ( i = 0 ; i < N ; i ++ ) { maxA = 0.0 ; imax = i ;               para ( k = i ; k < N ; k ++ ) si (( absA = fabs ( A [ k ][ i ])) > maxA ) { maxA = absA ; imax = k ; }                       si ( maxA < Tol ) devuelve 0 ; //falla, la matriz está degenerada       si ( imax != i ) { //pivotando P j = P [ i ]; P [ i ] = P [ imax ]; P [ imax ] = j ;               //pivotando filas de A ptr = A [ i ]; A [ i ] = A [ imax ]; A [ imax ] = ptr ;          //contando pivotes comenzando desde N (para determinante) P [ N ] ++ ; }   para ( j = i + 1 ; j < N ; j ++ ) { A [ j ][ i ] /= A [ i ][ i ];              para ( k = i + 1 ; k < N ; k ++ ) A [ j ][ k ] -= A [ j ][ i ] * A [ i ][ k ]; } }                 devuelve 1 ; //descomposición realizada }  /* ENTRADA: A,P rellenados en LUPDecompose; b - vector derecho; N - dimensión * SALIDA: x - vector solución de A*x=b */ void LUPSolve ( double ** A , int * P , double * b , int N , double * x ) {            para ( int i = 0 ; i < N ; i ++ ) { x [ i ] = b [ P [ i ]];             para ( int k = 0 ; k < i ; k ++ ) x [ i ] -= A [ i ][ k ] * x [ k ]; }               para ( int i = N - 1 ; i >= 0 ; i -- ) { para ( int k = i + 1 ; k < N ; k ++ ) x [ i ] -= A [ i ][ k ] * x [ k ];                            x [ i ] /= A [ i ][ i ]; } }   / * ENTRADA: A,P rellenados en LUPDecompose; dimensión N * SALIDA: IA es la inversa de la matriz inicial */ void LUPInvert ( double ** A , int * P , int N , double ** IA ) { for ( intj = 0 ; j < N ; j ++ ) { for ( inti = 0 ; i < N ; i ++ ) { IA [ i ] [ j ] = P [ i ] == j ? 1.0 : 0.0 ;                                        para ( int k = 0 ; k < i ; k ++ ) IA [ i ][ j ] -= A [ i ][ k ] * IA [ k ][ j ]; }               para ( int i = N - 1 ; i >= 0 ; i -- ) { para ( int k = i + 1 ; k < N ; k ++ ) IA [ i ][ j ] -= A [ i ][ k ] * IA [ k ][ j ];                            IA [ i ][ j ] /= A [ i ][ i ]; } } }    /* ENTRADA: A,P rellenados en LUPDecompose; dimensión N. * SALIDA: La función devuelve el determinante de la matriz inicial */ double LUPDeterminant ( double ** A , int * P , int N ) {        doble det = A [ 0 ][ 0 ];    para ( int i = 1 ; i < N ; i ++ ) det *= A [ i ][ i ];            devolver ( P [ N ] -N ) % 2 == 0 ?​ det : - det ; }           

Ejemplo de código C#

clase pública SystemOfLinearEquations { público double [] SolveUsingLU ( double [,] matriz , double [] parteDerecha , int n ) { // descomposición de la matriz double [,] lu = new double [ n , n ]; double suma = 0 ; para ( int i = 0 ; i < n ; i ++ ) { para ( int j = i ; j < n ; j ++ ) { suma = 0 ; para ( int k = 0 ; k < i ; k ++ ) suma += lu [ i , k ] * lu [ k , j ]; lu [ i , j ] = matriz [ i , j ] - suma ; } para ( int j = i + 1 ; j < n ; j ++ ) { suma = 0 ; para ( int k = 0 ; k < i ; k ++ ) suma += lu [ j , k ] * lu [ k , i ]; lu [ j , i ] = ( 1 / lu [ i , i ]) * ( matriz [ j , i ] - suma ); } }                                                                                                                   // lu = L+UI // encuentra la solución de Ly = b double [] y = new double [ n ]; for ( int i = 0 ; i < n ; i ++ ) { suma = 0 ; for ( int k = 0 ; k < i ; k ++ ) suma += lu [ i , k ] * y [ k ]; y [ i ] = rightPart [ i ] - suma ; } // encuentra la solución de Ux = y double [] x = new double [ n ]; for ( int i = n - 1 ; i >= 0 ; i -- ) { suma = 0 ; for ( int k = i + 1 ; k < n ; k ++ ) suma += lu [ i , k ] * x [ k ]; x [ i ] = ( 1 / lu [ i , i ]) * ( y [ i ] - suma ); } return x ; } }                                                                                            

Ejemplo de código MATLAB

función  LU = LUDecompDoolittle ( A ) n = length ( A ); LU = A ; % descomposición de la matriz, método de Doolittle para i = 1 : 1 : n para j = 1 :( i - 1 ) LU ( i , j ) = ( LU ( i , j ) - LU ( i , 1 :( j - 1 )) * LU ( 1 :( j - 1 ), j )) / LU ( j , j ); fin j = i : n ; LU ( i , j ) = LU ( i , j ) - LU ( i , 1 :( i - 1 )) * LU ( 1 :( i - 1 ), j ); fin %LU = L+UI fin                                             función  x = SolveLinearSystem ( LU, B ) n = longitud ( LU ); y = ceros ( tamaño ( B )); % encuentra la solución de Ly = B para i = 1 : n y ( i ,:) = B ( i ,:) - LU ( i , 1 : i ) * y ( 1 : i ,:); fin % encuentra la solución de Ux = y x = ceros ( tamaño ( B )); para i = n :( - 1 ): 1 x ( i ,:) = ( y ( i ,:) - LU ( i ,( i + 1 ): n ) * x (( i + 1 ): n ,:)) / LU ( i , i ); fin fin                                       A = [ 4 3 3 ; 6 3 3 ; 3 4 3 ] LU = LUDecompDoolittle ( A ) B = [ 1 2 3 ; 4 5 6 ; 7 8 9 ; 10 11 12 ] ' x = SolveLinearSystem ( LU , B ) A * x                                  

Véase también

• Descomposición del bloque LU
• Descomposición de Bruhat
• Descomposición de Cholesky
• Descomposición de la matriz de Crout
• Factorización LU incompleta
• Reducción de LU
• Descomposición matricial
• Descomposición QR

Notas

1. Schwarzenberg-Czerny, A. (1995). "Sobre la factorización matricial y la solución eficiente por mínimos cuadrados". Serie de suplementos de astronomía y astrofísica . 110 : 405. Código Bibliográfico :1995A&AS..110..405S.
2. Dwyer (1951).
3. Okunev y Johnson (1997), Corolario 3.
4. Trefethen y Bau (1997), pág. 166.
5. Trefethen y Bau (1997), pág. 161.
6. Lay, Lay y McDonald (2021), pág. 133, 2.5: Factorizaciones matriciales.
7. Rigotti (2001), Principal menor líder.
8. Horn & Johnson (1985), Corolario 3.5.5
9. Horn y Johnson (1985), Teorema 3.5.2.
10. Nhiayi, Ly; Phan-Yamada, Tuyetdong (2021). "Examinando la posible descomposición de LU". Revista GeoGebra de América del Norte . 9 (1).
11. Okunev y Johnson (1997).
12. Jefe de familia (1975).
13. Golub y Van Loan (1996), págs. 112, 119.
14. Cormen et al. (2009), p. 819, 28.1: Resolución de sistemas de ecuaciones lineales.
15. Shabat, Gil; Shmueli, Yaniv; Aizenbud, Yariv; Averbuch, Amir (2016). "Descomposición LU aleatoria". Análisis armónico computacional y aplicado . 44 (2): 246–272. arXiv : 1310.7202 . doi :10.1016/j.acha.2016.04.006. S2CID  1900701.
16. Bunch y Hopcroft (1974).
17. Trefethen y Bau (1997), pág. 152.
18. Golub y Van Loan (1996), pág. 121.

Referencias

• Bunch, James R.; Hopcroft, John (1974), "Factorización triangular e inversión mediante multiplicación rápida de matrices", Mathematics of Computation , 28 (125): 231–236, doi : 10.2307/2005828 , hdl : 1813/6003 , ISSN  0025-5718, JSTOR  2005828.
• Cormen, Thomas H. ; Leiserson, Charles E. ; Rivest, Ronald L. ; Stein, Clifford (2009), Introducción a los algoritmos (3.ª ed.), MIT Press y McGraw-Hill, ISBN 978-0-262-03293-3.
• Dwyer, Paul S. (1951), Cálculos lineales , Nueva York: Wiley.
• Golub, Gene H. ; Van Loan, Charles F. (1996), Cálculos matriciales (3.ª ed.), Baltimore: Johns Hopkins, ISBN 978-0-8018-5414-9.
• Horn, Roger A.; Johnson, Charles R. (1985), Análisis de matrices , Cambridge University Press, ISBN 978-0-521-38632-6. Véase la Sección 3.5. N  − 1
• Householder, Alston S. (1975), La teoría de matrices en el análisis numérico , Nueva York: Dover Publications , MR  0378371.
• Lay, David C.; Lay, Steven R.; McDonald, Judi J. (2021), Álgebra lineal y sus aplicaciones (sexta edición), Pearson, ISBN 978-0-13-585125-8.
• Okunev, Pavel; Johnson, Charles R. (1997), Condiciones necesarias y suficientes para la existencia de la factorización LU de una matriz arbitraria , arXiv : math.NA/0506382.
• Poole, David (2006), Álgebra lineal: una introducción moderna (2.ª ed.), Canadá: Thomson Brooks/Cole, ISBN 978-0-534-99845-5.
• Press, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007), "Sección 2.3", Recetas numéricas: el arte de la computación científica (3.ª ed.), Nueva York: Cambridge University Press, ISBN 978-0-521-88068-8.
• Trefethen, Lloyd N. ; Bau, David (1997), Álgebra lineal numérica , Filadelfia: Sociedad de Matemáticas Industriales y Aplicadas, ISBN 978-0-89871-361-9.
• Rigotti, Luca (2001), ECON 2001 - Introducción a los métodos matemáticos, Conferencia 8

Enlaces externos

Referencias

• Descomposición LU en MathWorld .
• Descomposición LU en Math-Linux .
• Descomposición LU en el Instituto de Métodos Numéricos Holísticos
• Factorización de matrices LU. Referencia MATLAB.

Código de computadora

• LAPACK es una colección de subrutinas FORTRAN para resolver problemas de álgebra lineal densos
• ALGLIB incluye un puerto parcial de LAPACK a C++, C#, Delphi, etc.
• Código C++, Prof. J. Loomis, Universidad de Dayton
• Código C, Biblioteca de fuentes de matemáticas
• Código de óxido
• LU en X10

Recursos en línea

• Aplicación web que resuelve de forma descriptiva sistemas de ecuaciones lineales con descomposición LU
• Calculadora de matrices con pasos, incluida la descomposición LU,
• Herramienta de descomposición LU, uni-bonn.de
• Descomposición LU por Ed Pegg, Jr. , The Wolfram Demonstrations Project , 2007.