Descomposición LU

Tipo de factorización matricial

En análisis numérico y álgebra lineal , la descomposición o factorización inferior-superior ( LU ) factoriza una matriz como el producto de una matriz triangular inferior y una matriz triangular superior (ver descomposición matricial ). En ocasiones, el producto 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 mediante 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] sólo 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 relación con problemas no simétricos... Banachiewicz... vio el punto... que el problema básico es en realidad uno de factorización matricial, o “descomposición”, como él llamó él." ^[2] A veces también se la conoce como descomposición LR (factores en matrices triangulares izquierda y derecha).

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 adecuadas, 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 por encima de la diagonal son cero, en la matriz triangular superior, todos los elementos 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, es posible que la factorización no se materialice. Por ejemplo, es fácil verificar (expandiendo la multiplicación de matrices ) que . Si , entonces al menos uno de y tiene que ser cero, lo que implica que L o U son singulares . Esto es imposible si A es no singular (invertible). Este es un problema de procedimiento. Se puede eliminar simplemente reordenando las filas de A de modo que el primer elemento de la matriz permutada sea distinto de cero. El mismo problema en pasos de factorización posteriores se puede eliminar 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 pivote parcial

Resulta que una permutación adecuada en filas (o columnas) es suficiente para la factorización LU. La factorización LU con pivote 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 por 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 de LUP sea una técnica útil en la práctica.

Factorización LU con pivote completo

Una factorización LU con pivote completo implica permutaciones de filas y 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 inferior-diagonal-superior (LDU)

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

${\displaystyle A=LDU,}$

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

Matrices rectangulares

Anteriormente requerimos que A sea una matriz cuadrada, pero todas estas descomposiciones también pueden generalizarse a matrices rectangulares. ^[6] En ese caso, L y D son matrices cuadradas que tienen el mismo número de filas que A , y U tiene exactamente las mismas dimensiones que A. Se debe interpretar que el triángulo superior 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 mediante 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á subdeterminado . En este caso, dos elementos cualesquiera distintos de cero de las matrices L y U son parámetros de la solución y se pueden establecer arbitrariamente en cualquier valor distinto de cero. Por lo tanto, para encontrar la descomposición LU única, es necesario imponer alguna restricción a las matrices L y U. Por ejemplo, podemos requerir convenientemente que la matriz triangular inferior L sea una matriz triangular unitaria (es decir, establecer todas las entradas de su diagonal principal en unidades). 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 unicidad

Matrices cuadradas

Cualquier 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 principales ^[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 menores principales principales son distintos de cero, aunque lo contrario 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 ) esté formada por unos. ${\textstyle L}$ ${\textstyle U}$

En general, cualquier matriz cuadrada podría tener una de las siguientes características: ${\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 cualquiera 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 son distintas de cero y linealmente independientes y al menos un principal menor principal es cero.

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

Matrices simétricas definidas positivas

Si A es una matriz definida positiva simétrica (o hermitiana , si A es compleja) , podemos arreglar las cosas de manera que U sea la transpuesta conjugada de L. Es decir, podemos escribir A como

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

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

Matrices generales

Para una matriz (no necesariamente invertible) sobre cualquier campo, 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 determinadas submatrices. El algoritmo de eliminación gaussiano 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, existe una fórmula cerrada (explícita) para los elementos de L , D y U en términos de proporciones de determinantes de ciertas submatrices de la matriz original A. ^[12] En particular, y para , es la relación entre la -ésima submatriz principal y la -ésima submatriz principal. El cálculo de los determinantes es 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)}$

Usando 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 pivote parcial añade sólo un término cuadrático; este no es el caso del giro completo. ^[13] ${\displaystyle {\tfrac {2}{3}}n^{3}}$

Explicación generalizada

Notación

Dada una matriz N × N , defina 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 han sido eliminados 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 y ayudarnos 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 esto, 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 matriz . Para nuestra matriz , podemos comenzar intercambiando filas para proporcionar las condiciones deseadas para la enésima columna. Por ejemplo, podríamos intercambiar filas para realizar un pivote parcial, o podríamos hacerlo para establecer el elemento de 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 todos los elementos siguientes en cero (donde está el elemento en la enésima columna de la diagonal principal). Denotaremos cada elemento a continuación como (dónde ). Para establecerlo en cero, configuramos 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 denominada , ingresando directamente los valores calculados previamente mediante la siguiente fórmula. ${\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 un pivote parcial y así intercambiar la primera y segunda fila para que nuestra matriz y la primera iteración de nuestra matriz respectivamente se conviertan en Una vez que hayamos intercambiado las filas, podemos eliminar los elementos debajo de la diagonal principal. en la primera columna actuando de manera que, una vez restadas estas filas, hemos derivado de la matriz. Debido a que estamos implementando un pivote parcial, intercambiamos la segunda y 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 tal que . Debido a que no existen elementos distintos de cero debajo de la diagonal principal en nuestra iteración actual después de esta resta de filas, esta resta de filas deriva nuestra matriz final (denotada como ) y nuestra matriz final: Después de cambiar 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 dónde está la matriz identidad N × N con su n -ésima columna reemplazada por el vector transpuesto . En otras palabras, el triángulo inferior matriz ${\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 usando la fórmula equivale 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 unos 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, multiplique y genere la matriz fusionada indicada 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 momento, es necesario intercambiar la n -ésima fila con otra fila debajo antes de continuar. Esta es la razón por la que en general parece una descomposición LU . ${\displaystyle a_{n,n}^{(n-1)}\neq 0}$ ${\displaystyle \ell _{i,n}}$ ${\displaystyle P^{-1}A=LU}$

Descomposición de LU Crout

Tenga en cuenta que la descomposición obtenida mediante este procedimiento es una descomposición de Doolittle : la diagonal principal de L se compone únicamente de 1 s. Si se procediera a eliminar elementos encima de la diagonal principal sumando múltiplos de las columnas (en lugar de eliminar elementos debajo de la diagonal sumando múltiplos de las filas ), obtendríamos 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 A dada es obtener una descomposición de Doolittle de la transpuesta de A. De hecho, si la descomposición LU se obtiene mediante el algoritmo presentado en esta sección, entonces, tomando y , tenemos que se trata de 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 recursividad

Cormen et al. ^[14] describen un algoritmo recursivo para la descomposición de 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 matriz identidad en caso contrario. Ahora vamos , si ; o de otro modo. 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 . Dejar . 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, de manera que con alta probabilidad , donde es una constante que depende de los parámetros del algoritmo y es la -ésima 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 de matriz dispersa

Se han desarrollado algoritmos especiales para factorizar matrices dispersas grandes . Estos algoritmos intentan encontrar factores dispersos L y U. Idealmente, el costo del cálculo está determinado por el número de entradas distintas de cero, en lugar del 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 mediante 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. Segundo, resolvemos la ecuación para x . ${\textstyle U\mathbf {x} =\mathbf {y} }$

En ambos casos estamos tratando con 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 gaussiano (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 las 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 gaussiano.

El costo de resolver un sistema de ecuaciones lineales es aproximadamente operaciones de punto flotante si la matriz tiene tamaño . Esto lo hace dos veces más rápido que los algoritmos basados ​​en la descomposición QR , que cuestan aproximadamente operaciones de punto flotante cuando se utilizan reflexiones de Householder . Por esta razón, normalmente se prefiere la descomposición LU. ^[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 generalmente se trata como 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 de n por p , de modo que estamos tratando de encontrar una matriz X (también una matriz de n por p ):

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

Podemos usar 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 seguiría que el resultado X debe ser el inverso de A. ^[18]

Calcular 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 pivote 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 haciendo que P sea 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 - número de tolerancia pequeño para detectar fallas cuando la matriz está casi degenerada * SALIDA: La matriz A se cambia, 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 , doble Tol , int * P ) {          int yo , 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 ; imáx = yo ;               for ( k = i ; k < N ; k ++ ) if (( absA = fabs ( A [ k ][ i ])) > maxA ) { maxA = absA ; imáx = k ; }                       si ( maxA < Tol ) devuelve 0 ; //fallo, la matriz está degenerada       if ( imax ! = i ) { //pivotando P j = P [ i ]; P [ i ] = P [ imax ]; P [ imax ] = j ;               //filas pivotantes de A ptr = A [ i ]; A [ i ] = A [ imax ]; A [ imax ] = ptr ;          //contando pivotes a partir de N (para determinante) P [ N ] ++ ; }   for ( 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 ]; } }                 devolver 1 ; //descomposición realizada }  /* ENTRADA: A,P completado LUPDecompose; b - vector rhs; 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 ]; }               for ( int i = N - 1 ; i >= 0 ; i -- ) { for ( int k = i + 1 ; k < N ; k ++ ) x [ i ] -= A [ i ][ k ] * x [ k ];                            x [ yo ] /= A [ yo ][ yo ]; } }   /* ENTRADA: A,P completado LUPDecompose; N - dimensión * SALIDA: IA es la inversa de la matriz inicial */ void LUPInvert ( double ** A , int * P , int N , double ** IA ) { for ( int j = 0 ; j < N ; j + + ) { for ( int i = 0 ; i < N ; i ++ ) { IA [ i ][ j ] = P [ i ] == j ? 1,0 : 0,0 ;                                        for ( int k = 0 ; k < i ; k ++ ) IA [ i ][ j ] -= A [ i ][ k ] * IA [ k ][ j ]; }               for ( int i = N - 1 ; i >= 0 ; i -- ) { for ( 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 completado LUPDecompose; norte - dimensió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 ];    for ( 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 doble [] SolveUsingLU ( doble [,] matriz , doble [] parte derecha , int n ) { // descomposición de la matriz doble [,] lu = nuevo doble [ n , n ]; doble 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 ]; para ( int i = 0 ; i < n ; i ++ ) { suma = 0 ; para ( int k = 0 ; k < i ; k ++ ) suma += lu [ i , k ] * y [ k ]; y [ i ] = partederecha [ i ] -suma ;​ } // encuentra la solución de Ux = y double [] x = new double [ n ]; for ( int i = n - 1 ; i >= 0 ; i - ) { suma = 0 ; para ( int k = i + 1 ; k < n ; k ++ ) suma += lu [ i , k ] * x [ k ]; x [ i ] = ( 1 / lu [ i , i ]) * ( y [ i ] - suma ); } devolver x ; } }                                                                                            

Ejemplo de códigoMATLAB

función  LU = LUDecompDoolittle ( A ) n = longitud ( A ); LU = A ; % de 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 ); final j = yo : norte ; LU ( i , j ) = LU ( i , j ) - LU ( i , 1 : (i - 1 )) * LU ( 1 : (i - 1 ), j ); final %LU = L+UI final                                             función  x = SolveLinearSystem ( LU, B ) n = longitud ( LU ); y = ceros ( tamaño ( B )); % encontrar solución de Ly = B para i = 1 : n y ( i ,:) = B ( i ,:) - LU ( i , 1 : i ) * y ( 1 : i ,:); final % encontrar 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 ( yo , yo ); 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                                  

Ver también

• Descomposición de LU en bloque
• 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 factorización matricial y solución eficiente de mínimos cuadrados". Serie de Suplementos de Astronomía y Astrofísica . 110 : 405. Código bibliográfico : 1995A y AS..110..405S.
2. Dwyer, Paul S. (1951). Cálculos lineales . Nueva York: Wiley.
3. Okunev & Johnson (1997), Corolario 3.
4. Trefethen y Bau (1997), pág. 166.
5. Trefethen y Bau (1997), pág. 161.
6. Lay, David C. (2016). Álgebra lineal y sus aplicaciones. Steven R. Lay, Judith McDonald (Quinta ed.). Harlow. pag. 142.ISBN​ 978-1-292-09223-2. OCLC  920463015.{{cite book}}: CS1 maint: location missing publisher (link)
7. Rigotti (2001), Principal menor principal
8. Horn & Johnson (1985), Corolario 3.5.5
9. Horn & Johnson (1985), Teorema 3.5.2
10. Nhiayi, Ly; Phan-Yamada, Tuyetdong (2021). "Examen de la posible descomposición de LU". Revista norteamericana GeoGebra . 9 (1).
11. Okunev y Johnson (1997)
12. Jefe de familia (1975)
13. Préstamo Golub y Van (1996), pág. 112, 119.
14. Cormen, Thomas H .; Leiserson, Charles E .; Rivest, Ronald L .; Stein, Clifford (2001). Introducción a los algoritmos . MIT Press y McGraw-Hill. ISBN 978-0-262-03293-3.
15. Shabat, Gil; Shmueli, Yaniv; Aizenbud, Yariv; Averbuch, Amir (2016). "Descomposición LU aleatoria". Análisis Armónico Aplicado y Computacional . 44 (2): 246–272. arXiv : 1310.7202 . doi :10.1016/j.acha.2016.04.006. S2CID  1900701.
16. Manojo y Hopcroft (1974)
17. Trefethen y Bau (1997), pág. 152.
18. Préstamo Golub y Van (1996), pág. 121

Referencias

• Manojo, James R.; Hopcroft, John (1974), "Factorización e inversión triangular mediante multiplicación rápida de matrices", Matemáticas de la Computación , 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 (2001), Introducción a los algoritmos , MIT Press y McGraw-Hill, ISBN 978-0-262-03293-3.
• Golub, Gene H .; Van Loan, Charles F. (1996), Computaciones matriciales (3.ª ed.), Baltimore: Johns Hopkins, ISBN 978-0-8018-5414-9.
• Cuerno, Roger A.; Johnson, Charles R. (1985), Análisis matricial , Cambridge University Press, ISBN 978-0-521-38632-6. Consulte la Sección 3.5. norte  -1
• Jefe de familia, Alston S. (1975), La teoría de las matrices en el análisis numérico , Nueva York: Dover Publications , MR  0378371.
• 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.
• Prensa, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007), "Sección 2.3", Recetas numéricas: el arte de la informática 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 matricial LU. Referencia de MATLAB.

Codigo de computadora

• LAPACK es una colección de subrutinas FORTRAN para resolver problemas de álgebra lineal densa
• 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 fuente de matemáticas
• código de óxido
• LU en X10

Recursos en línea

• WebApp resolución descriptiva de sistemas de ecuaciones lineales con Descomposición LU
• Calculadora matricial 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.