descomposición LU

En análisis numérico y álgebra lineal , la descomposición o factorización LU ( de matriz inferior a matriz superior ) factoriza una matriz como el producto de una matriz triangular inferior y una matriz triangular superior (véase multiplicación de matrices y descomposición de matrices ). El producto a veces incluye también una matriz de permutación . La descomposición LU puede considerarse como la forma matricial de la eliminación gaussiana . Los ordenadores suelen resolver sistemas de ecuaciones lineales cuadradas mediante la descomposición LU, y también es un paso clave al invertir una matriz o calcular su determinante . A veces también se la denomina descomposición LR (factoriza en matrices triangulares izquierdas y derechas). El algoritmo de descomposición LU para matrices generales fue introducido por el astrónomo polaco Tadeusz Banachiewicz en 1938. ^[¹^]

Definiciones

Sea $A$ una matriz cuadrada. Una factorización LU se refiere a la expresión de $A$ en el producto de dos factores: una matriz triangular inferior $L$ y una matriz triangular superior $U,$ de modo que $A = LU$ . A veces, la factorización es imposible sin una reordenación previa de $A$ para evitar la división por cero o el crecimiento incontrolado de errores de redondeo. Por lo tanto, la expresión alternativa es $PAQ = LU$ , donde en notación formal los factores de matriz de permutación $P$ y $Q$ indican la permutación de filas (o columnas) de $A.$ En teoría, $P$ (o $Q$ ) se obtiene mediante permutaciones de filas (o columnas) de la matriz identidad ; en la práctica , las permutaciones correspondientes se aplican directamente a las filas (o columnas) de $A.$

Matriz $A$ de lado $n$ tiene $n^{2}$ Los coeficientes de dos matrices triangulares combinadas contienen $n (n +1)$ coeficientes, y por lo tanto, $los n$ coeficientes de las matrices $LU$ no son independientes. La convención es establecer $L$ como unitarangular, es decir, con todos $los n$ elementos de la diagonal principal iguales a uno. Sin embargo, establecer $U$ como unitarangular reduce el procedimiento al mismo después de la transposición del producto matricial (véanse las propiedades de la transposición matricial): $B=A^{\textsf {T}}=(LU)^{\textsf {T}}=U^{\textsf {T}}L^{\textsf {T}}.$ Tras la transposición, $U T$ es el triángulo inferior, mientras que $L T$ es el factor unitarigular superior de $B.$ Esto demuestra también que las operaciones sobre filas (por ejemplo, el pivoteo) son equivalentes a las realizadas sobre columnas de una matriz transpuesta, y que, en general, la elección del algoritmo de filas o columnas no ofrece ninguna ventaja.

En la matriz triangular inferior, todos los elementos por encima de la diagonal principal 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í: ${\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 o permutaciones adecuadas en la matriz, la factorización puede no materializarse. Por ejemplo, es fácil verificar (expandiendo la multiplicación de matrices ) que ${\textstyle a_{11}=\ell _{11}u_{11}}$ . Si ${\textstyle a_{11}=0}$ , entonces al menos uno de ${\textstyle \ell _{11}}$ y ${\textstyle u_{11}}$ tiene que ser cero, lo que implica que $L$ o $U$ es singular . Esto es imposible si $A$ es no singular (invertible). En términos de operaciones, la anulación/eliminación de los elementos restantes de la primera columna de $A$ implica la división de ${\textstyle a_{21},a_{31}}$ con ${\textstyle a_{11}}$ , imposible si es 0. 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 los pasos de factorización subsiguientes se puede eliminar de la misma manera. Para la estabilidad numérica frente a errores de redondeo/división por números pequeños es importante seleccionar ${\textstyle a_{11}}$ de gran valor absoluto (cf. pivotar).

LU mediante recursión

El ejemplo anterior de matrices de 3 × 3 demuestra que el producto matricial de la fila superior y las columnas más a la izquierda de las matrices involucradas juega un papel especial para que $LU$ tenga éxito. Marquemos versiones consecutivas de matrices con $(0),\;(1),\dots$ y luego escribamos el producto matricial $A\equiv A^{(0)}=L^{(0)}U^{(0)}$ de tal manera que estas filas y columnas queden separadas del resto. Para ello, utilizaremos la notación de matriz de bloques , de modo que, por ejemplo, $a\equiv a_{11}$ es un número ordinario, ${\bf {w}}^{\textsf {T}}\equiv (a_{12},a_{13})^{\textsf {T}}$ es un vector fila y ${\bf {v}}=(a_{21},a_{31})$ es un vector columna y $A'$ es una submatriz de la matriz $A^{(0)}$ sin la fila superior y la columna más a la izquierda. Entonces podemos reemplazar $A^{(0)}=L^{(0)}U^{(0)}$ con un producto de bloques de matrices . Es decir, resulta que se pueden multiplicar bloques de matrices como si fueran números ordinarios, es decir, fila por columna, excepto que ahora sus componentes son submatrices, a veces reducidas a escalares o vectores. Por lo tanto, $u{\bf {l}}$ denota un vector obtenido de ${\bf {l}}$ después de multiplicar cada componente por un número $u$ , ${\bf {lu}}^{\textsf {T}}$ es un producto exterior de vectores ${\bf {l,u}}$ , es decir, una matriz cuya primera columna es $u_{12}{\bf {l}}$ , el siguiente es $u_{13}{\bf {l}}$ y así sucesivamente para todos los componentes de ${\bf {u}}$ y $L^{(1)}U^{(1)}$ es un producto de submatrices de $L^{(0)},\;U^{(0)}$ ${\begin{aligned}\left({\begin{array}{c|c}a&{\bf {w}}^{\textsf {T}}\\\hline \\[-0.5em]{\bf {v}}&\quad A'\quad \\[-0.5em]\\\end{array}}\right)&=\left({\begin{array}{c|c}{\rm {1}}&{\bf {0}}^{\textsf {T}}\\\hline \\[-0.5em]{\bf {l}}&\quad L^{(1)}\quad \\[-0.5em]\\\end{array}}\right)\;\left({\begin{array}{c|c}u&{\bf {u}}^{\textsf {T}}\\\hline \\[-0.5em]{\bf {0}}&\quad U^{(1)}\\[-0.5em]\\\end{array}}\right)\\&=\left({\begin{array}{c|c}u&{\bf {u}}^{\textsf {T}}\\\hline \\[-0.5em]u{\bf {l}}&\quad {\bf {lu}}^{\textsf {T}}+L^{(1)}U^{(1)}\\[-0.5em]\\\end{array}}\right)\end{aligned}}$

De la igualdad de la primera y la última matriz se deduce la final. $u=a$ , ${\bf {u=w}}$ , ${\bf {l}}={(1/a)}{\bf {v}}$ mientras que matriz $A'$ se actualiza/reemplaza con $A^{(1)}\equiv L^{(1)}U^{(1)}=$ $A'-{\bf {lu}}^{\textsf {T}}$ . Ahora viene la observación crucial: nada nos impide tratar $A^{(1)}$ de la misma manera que lo hicimos con $A^{(0)}$ , repetidamente. Si la dimensión de $A$ es $n$ × $n$ , después de $n - 1$ pasos de este tipo todas las columnas ${\bf {v}}$ forman la parte subdiagonal de la matriz triangular $L$ y todos los pivotes $a$ combinado con filas ${\bf {w}}^{\textsf {T}}$ formar matriz triangular superior $U$ , según se requiera. En el ejemplo anterior $n = 3$ , por lo que solo dos pasos son suficientes.

El procedimiento anterior demuestra que en ningún paso el elemento pivote diagonal superior $a$ de submatrices consecutivas pueden ser cero. Para evitarlo, se pueden intercambiar columnas o filas de modo que $a$ se vuelve distinto de cero. Este procedimiento que implica permutación se llama LUP , descomposición con pivoteo.

La permutación de columnas corresponde al producto de la matriz $AQ^{(0)}$ dónde $Q^{(0)}$ es una matriz de permutación, es decir, la matriz identidad. $I$ después de la misma permutación de columna. Después de todos los pasos, dicha descomposición LUP se aplica a $AQ^{(0)}\cdots Q^{(n-1)}\equiv AQ=LU$ El esquema de cálculo actual y otros similares en Cormen et al. ^{[ 2 ]} son ejemplos de algoritmos de recurrencia . Demuestran dos propiedades generales de la factorización LU:

la necesidad de adaptarse en cada paso; y
Los valores finales de las matrices $L$ y $U$ se obtienen gradualmente, una fila o una columna por paso.

Los algoritmos de recurrencia no son excesivamente costosos en términos de operaciones algebraicas, pero presentan la desventaja práctica de tener que actualizar y almacenar la mayoría de los elementos de $A$ en cada paso. Se verá que, al reordenar los cálculos, es posible prescindir del almacenamiento de valores intermedios.

Factorización LU con pivoteo parcial

Resulta que una permutación adecuada de filas (o columnas) para seleccionar el pivote máximo absoluto de columna (o fila) $a 11$ es suficiente para una factorización LU numéricamente estable, excepto en casos patológicos conocidos. Se denomina « factorización LU con pivoteo parcial » (LUP). $PA=LU,\quad (AQ=LU),$ donde $L$ y $U$ son nuevamente matrices triangulares inferior y superior, y $P$ y $Q son$ matrices de permutación correspondientes que, al multiplicarse por la izquierda y por la derecha respectivamente por $A$ , reordenan las filas y columnas de $A.$ Resulta que todas las matrices cuadradas pueden factorizarse 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.

Una variante denominada " pivote de torre " implica, en cada paso, la búsqueda del elemento máximo, siguiendo el movimiento de una torre en un tablero de ajedrez: columna, fila, columna nuevamente, y así sucesivamente, hasta alcanzar un pivote que sea máximo tanto en su fila como en su columna. Se puede demostrar que, para matrices grandes de elementos aleatorios, el costo de las operaciones en cada paso es similar al del pivote parcial, proporcional a la longitud del lado de la matriz, a diferencia del pivote completo, que es proporcional al cuadrado de dicha longitud.

Factorización LU con pivoteo completo

Una ' factorización LU con pivoteo completo ' implica permutaciones tanto de filas como de columnas para encontrar el elemento máximo absoluto en toda la submatriz: $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 ' (DDI) es una descomposición de la forma $A=LDU,$ donde $D$ es una matriz diagonal , y $L$ y $U$ son matrices unitarangulares , lo que significa que todas las entradas en las diagonales de $L$ y $U$ son uno.

Matrices rectangulares

Arriba requerimos que $A$ sea una matriz cuadrada, pero estas descomposiciones también se pueden generalizar 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.$ 'Triangular superior' debe interpretarse como tener solo entradas cero 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 × 2 : ${\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 $\left\{{\begin{alignedat}{4}\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{alignedat}}\right.$

Este sistema de ecuaciones es subdeterminado . En este caso, cualesquiera dos elementos no nulos de las matrices $L$ y $U$ son parámetros de la solución y pueden asignarse arbitrariamente a 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 exigir convenientemente que la matriz triangular inferior $L$ sea una matriz triangular unitaria, de modo que todos los elementos de su diagonal principal sean iguales a uno. Entonces, el sistema de ecuaciones tiene la siguiente solución: ${\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 ${\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

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

Si una matriz cuadrada e invertible tiene una factorización LDU (con todos los elementos 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 $L$ o $U$ esté formada por unos.

En general, cualquier matriz cuadrada $A n \times n$ podría tener una de las siguientes características:

una factorización LU única (como se mencionó anteriormente);
infinitas factorizaciones LU si alguna de las primeras $(n -1)$ columnas es linealmente dependiente;
No hay factorización LU si las primeras $(n - 1)$ columnas son linealmente independientes y al menos un menor principal principal es cero.

En el caso 3, se puede aproximar una factorización LU cambiando una entrada diagonal $a ij$ por $a ij \pm ε$ para evitar un menor principal principal cero. ^{[ 10 ]}

Matrices simétricas definidas positivas

Si $A$ es una matriz simétrica (o hermitiana , si $A$ es compleja) definida positiva , podemos ordenar las cosas de manera que $U$ sea la transpuesta conjugada de $L.$ Es decir, podemos escribir $A$ como $A=LL^{*}\,.$

Esta descomposición se denomina descomposición de Cholesky . Si $A$ 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.

Matrices generales

Para una matriz (no necesariamente invertible) sobre cualquier cuerpo, se conocen las condiciones exactas necesarias y suficientes para que tenga una factorización LU. Estas 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 ú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, $D 1 = A 1,1$ , y para $i = 2, ..., n$ , $D i$ es la razón de la $i$ -ésima submatriz principal a la $(i - 1)$ -ésima submatriz principal. El cálculo de los determinantes es computacionalmente costoso , por lo que esta fórmula explícita no se usa en la práctica.

Utilizando la eliminación gaussiana

El siguiente algoritmo es esencialmente una forma modificada de eliminación gaussiana . El cálculo $de$ una descomposición LU mediante este algoritmo requiere $2/3 n³ operaciones de$ punto flotante, ignorando los términos de orden inferior. El pivoteo parcial $añade$ $solo$ un término cuadrático; esto no ocurre con ^el pivoteo completo. $[$ ^{13 ]}

Explicación generalizada

Notación

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

A continuación se muestra una matriz para observar y así ayudarnos a recordar la notación (donde cada $*$ representa cualquier número real en la matriz):

$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 $A$ mediante operaciones de fila hasta que se convierte en la matriz $U,$ en la que todos los elementos por debajo de la diagonal principal son iguales a cero. Durante este proceso, crearemos simultáneamente dos matrices separadas, $P$ y $L$ , de modo que $PA = LU$ .

Definimos la matriz de permutación final $P$ como la matriz identidad, cuyas filas se intercambian en el mismo orden que la matriz $A$ al transformarse en la matriz $U.$ Para nuestra matriz $A (n -1)$ , podemos comenzar intercambiando filas para proporcionar las condiciones deseadas para la columna $n$ . Por ejemplo, podríamos intercambiar filas para realizar un pivoteo parcial, o podríamos hacerlo para establecer el elemento pivote $a n,n$ en la diagonal principal a un número distinto de cero para poder completar la eliminación gaussiana.

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

También podemos crear la matriz triangular inferior denotada como $L$ , introduciendo directamente los valores previamente calculados de $ℓ i,n$ mediante la fórmula que aparece a continuación.

$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 se nos da la matriz $A={\begin{pmatrix}0&5&{\frac {22}{3}}\\4&2&1\\2&7&9\\\end{pmatrix}},$ Optaremos por implementar el pivoteo parcial y, por lo tanto, intercambiaremos la primera y la segunda fila de manera que nuestra matriz $A$ y la primera iteración de nuestra matriz $P, respectivamente, se conviertan en...$ $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}}.$ Una vez que hayamos intercambiado las filas, podemos eliminar los elementos que se encuentran debajo de la diagonal principal en la primera columna realizando ${\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}}$ de tal manera que, ${\begin{alignedat}{0}\ell _{2,1}={\frac {0}{4}}=0\\\ell _{3,1}={\frac {2}{4}}=0.5\end{alignedat}}$ Una vez restadas estas filas, hemos obtenido de $A (1)$ la matriz $A^{(1)}={\begin{pmatrix}4&2&1\\0&5&{\frac {22}{3}}\\0&6&8.5\\\end{pmatrix}}.$

Debido a que estamos implementando un pivoteo parcial, intercambiamos la segunda y tercera fila de nuestra matriz derivada y la versión actual de nuestra matriz $P$ respectivamente para obtener $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}}.$ Ahora, eliminamos los elementos debajo de la diagonal principal en la segunda columna realizando $la fila 3 = fila 3 - (ℓ 3,2)\cdot fila 2$ tal que $ℓ 3,2 = ⁠ 5 / 6 ⁠$ . Debido a que no existen elementos distintos de cero debajo de la diagonal principal en nuestra iteración actual de $A$ después de esta resta de filas, esta resta de filas deriva nuestra matriz $A$ final (denotada como $U ) y la matriz$ $P$ final : $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}}.$ Tras intercambiar también las filas correspondientes, obtenemos nuestra matriz $L final:$ $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}}$

Ahora bien, estas matrices tienen una relación tal que $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 $n$ estableciendo $A^{(n)}:=L_{n}^{-1}A^{(n-1)},$ dónde $L_{n}^{-1}$ es la matriz identidad $N$ × $N$ con su $n$ -ésima columna reemplazada por el vector transpuesto $(0 \dots 01 - ℓ n +1, n \dots - ℓ N, n) T$ .

En otras palabras, la matriz triangular inferior $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 $N - 1$ columnas utilizando el $A^{(n)}:=L_{n}^{-1}A^{(n-1)}$ La fórmula es equivalente a encontrar la descomposición. $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)}.$ Denota $L = L 1 \dots L N -1$ de modo que $A = LA (N -1) = LU$ .

Ahora calculemos la secuencia de $L 1 \dots L N -1$ . Sabemos que $L i$ tiene la siguiente fórmula: $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 posición en el producto de las dos matrices. Por ejemplo:

$\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 $L 1$ juntos y generamos la matriz fusionada denotada como $L$ (como se mencionó anteriormente). Usando la matriz $L$ , obtenemos $A = LU$ .

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

Descomposición de LU Banachiewicz

Aunque el algoritmo de descomposición LU de Banachiewicz (1938) precedió a la llegada de las computadoras electrónicas programadas, estaba listo para su implementación directa en código, ya que el intercambio de índices, la transposición y la multiplicación columna por columna siguen siendo capacidades nativas de la mayoría de los lenguajes de programación y son manejadas únicamente por los compiladores con poco retraso en la ejecución real. La peculiar notación matricial utilizada por Banachiewicz le permitió multiplicar matrices columna por columna, una característica conveniente para cálculos mecánicos, ya que podía revelar factores consecutivos deslizando una regla a las siguientes filas de las matrices. Sin embargo, para los lectores humanos, sus ecuaciones se transforman mejor a la notación matricial estándar. Para obtener, a partir de una matriz completa $A,$ los cálculos de las matrices triangulares $U$ y $L$ , comience copiando la fila superior y la columna más a la izquierda de $A$ respectivamente en las posiciones correspondientes de las matrices $U$ y $L.$ Los elementos diagonales unitarios conocidos de $L$ no se almacenan ni se utilizan durante todo el proceso. Los siguientes cálculos continúan para las filas y columnas subsiguientes hasta la esquina inferior derecha de $A.$

La figura ilustra los cálculos para la tercera fila y columna, suponiendo que las etapas anteriores ya se completaron. Las matrices involucradas se nombran sobre los cuadrados que marcan su contenido. Los productos y restas de matrices se aplican solo a los elementos dentro de los recuadros gruesos. Los recuadros finos rellenos de verde indican valores ya conocidos, de etapas anteriores. Los recuadros azules indican los lugares en las matrices $U$ y $L$ donde se almacenan los resultados. Tenga en cuenta que en cada etapa, los elementos resultantes de $L$ deben dividirse por el elemento pivote correspondiente en la diagonal principal de $U.$ Esto también se aplica a la columna más a la izquierda de $L.$

Cabe destacar que, tras completar la tercera etapa, los elementos de la matriz $A$ ya no se utilizan, ni tampoco los de las etapas anteriores. Esto permite reemplazar dichos elementos con los valores resultantes de $U$ y $L$ , es decir, ejecutar la descomposición LU in situ , de modo que toda la matriz $A$ se reemplaza con $U$ y $L,$ excepto la diagonal unitaria de $L.$ El algoritmo LU de Banachiewicz es idóneo para el pivoteo parcial, ya que selecciona el pivote máximo absoluto de la fila recién calculada de $U$ y, posteriormente, intercambia sus columnas para que coincida con la diagonal principal. Se pueden obtener más detalles examinando el código Fortran90 adjunto.

Todos los algoritmos LU de pivote parcial cuestan aproximadamente la misma cantidad, de orden ${\textstyle O\left({2 \over 3}n^{3}\right)}$ operaciones, donde $n$ es el número de filas o columnas de $A.$

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 unos. Si procediéramos eliminando los elementos situados por encima de la diagonal principal mediante la suma de múltiplos de las columnas (en lugar de eliminar los elementos situados por debajo de la diagonal mediante la suma de múltiplos de las filas ), obtendríamos una descomposición de Crout , donde la diagonal principal de $U$ está compuesta por unos.

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 $A T = L 0 U 0$ es la descomposición LU obtenida a través del algoritmo presentado en esta sección, entonces al tomar $L = U T 0$ y $U = L T 0$ , tenemos que $A = LU$ es una descomposición de Crout.

Algoritmo aleatorio

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

Complejidad teórica

Si dos matrices de orden $n$ se pueden multiplicar en tiempo $M (n)$ , donde $M (n) \geq n a$ para algún $a > 2$ , entonces se puede calcular una descomposición LU en tiempo $O(M (n))$ . ^{[ 15 ]} Esto significa, por ejemplo, que existe un algoritmo $O(n²³⁷⁶)$ basado en el algoritmo de Coppersmith-Winograd . Véase también el artículo sobre algoritmos rápidos de multiplicación de matrices $para$ más detalles.

Descomposición de matrices dispersas

Se han desarrollado algoritmos especiales para factorizar matrices dispersas de gran tamaño . Estos algoritmos buscan los factores dispersos $L$ y $U.$ Idealmente, el costo computacional viene determinado por el número de entradas no nulas, 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 de cero a un valor distinto de cero durante la ejecución de un algoritmo).

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

Aplicaciones

Resolución de ecuaciones lineales

Dado un sistema de ecuaciones lineales en forma matricial $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 $PA = LU$ , entonces $LU x = P b$ .

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

Primero, resolvemos la ecuación $L y = P b$ para $y$ .
Segundo, resolvemos la ecuación $U x = y$ para $x$ .

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

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

$El$ coste de resolver un sistema de ecuaciones lineales es aproximadamente $2/3$ $n³ operaciones$ $de$ punto flotante si la matriz $A$ tiene tamaño $n .$ $Esto$ lo hace dos veces más rápido que los algoritmos basados en $la$ descomposición QR , que cuestan alrededor de $4/3$ $n³$ $operaciones$ $de$ $punto flotante$ cuando se utilizan reflexiones de Householder $. Por esta razón, generalmente$ $se$ $prefiere la descomposición LU.$ ^[¹⁶^]

Invertir una matriz

Al resolver sistemas de ecuaciones, $b$ se suele tratar 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 $n$ × $p$ , por lo que estamos tratando de encontrar una matriz $X$ (también una matriz $n$ × $p$ ): $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$ , $I n$ . De ello se deduce que el resultado $X$ debe ser la inversa de $A.$ ^[¹⁷ ]

Calcular el determinante

Dada la descomposición LUP $A = P -1 LU$ de una matriz cuadrada $A$ , el determinante de $A$ se puede calcular directamente como ${\begin{aligned}\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).\end{aligned}}$

La segunda ecuación se deduce 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, $det(A)$ también es igual al lado derecho de la ecuación anterior, si hacemos $S$ como el número total de intercambios de filas y columnas.

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

Historia

La descomposición LU está relacionada con la eliminación de sistemas de ecuaciones lineales, como describió Ralston. ^{[ 18 ]} La solución de $N$ ecuaciones lineales con $N$ incógnitas mediante eliminación ya era conocida por los antiguos chinos. ^{[ 19 ]} Antes de Gauss, muchos matemáticos en Eurasia la aplicaban y perfeccionaban, pero como el método se relegó al ámbito escolar, pocos dejaron descripciones detalladas. Por lo tanto, el nombre de eliminación gaussiana es solo una abreviatura conveniente de una historia compleja.

El astrónomo polaco Tadeusz Banachiewicz introdujo la descomposición LU en 1938. ^{[ 20 ]} Sobre Banachiewicz, Paul Dwyer afirmó: ^{[ 21 ]}

Parece ser que Gauss y Doolittle aplicaron el método [de eliminación] solo a ecuaciones simétricas. Autores más recientes, como Aitken, Banachiewicz, Dwyer y Crout , han hecho hincapié en el uso del método, o variaciones del mismo, en relación con problemas no simétricos . Banachiewicz comprendió que el problema fundamental es, en realidad, de factorización matricial, o "descomposición", como él la denominó.

— Paul Dwyer, Cálculos lineales (1951)

Banachiewicz ^{[ 20 ]} fue el primero en considerar la eliminación en términos de matrices y de esta manera formuló la descomposición LU, como lo demuestra su ilustración gráfica. Sus cálculos siguen los de las matrices ordinarias, pero la notación difiere en que prefirió escribir un factor transpuesto, para poder multiplicarlos mecánicamente columna por columna, deslizando una regla por filas consecutivas de ambos (usando un aritmómetro ). Combinado con el orden intercambiado de los índices, sus fórmulas en notación moderna se leen ${\begin{aligned}{\mathbf {x} }\cdot IA'={\mathbf {0} }&\rightarrow A'{\mathbf {x} =0}\equiv (A|{\mathbf {l} }){\mathbf {x} },\\A=G\cdot H&\rightarrow A^{T}=G^{T}H,\end{aligned}}$

donde $IA \to A T$ ; $x \equiv [x 1, ... , x n, -1]$ ; $A'$ se refiere a $A$ extendida con la última columna; y el último componente de $x$ es $-1$ . Las fórmulas matriciales para calcular filas y columnas de factores LU por recursión se dan en la parte restante del artículo de Banachiewicz como Eq. (2.3) y (2.4) . Este artículo de Banachiewicz contiene tanto la derivación de factores $LU$ como $R T R$ de matrices no simétricas y simétricas respectivamente. A veces se confunden ya que publicaciones posteriores tienden a vincular su nombre únicamente con el redescubrimiento de la descomposición de Cholesky. El propio Banachiewicz puede ser excusado de inacción ya que al año siguiente sufrió persecución por parte de los ocupantes, pasando tres meses en el campo de concentración de Sachsenhausen , de donde, al ser liberado, llevó consigo desde un tren a su colaborador y compañero de prisión Antoni Wilk, quien murió de agotamiento una semana después.

Ejemplos de código

Ejemplo de código Fortran90

Módulo mlu Implícito Ninguno Entero , Parámetro :: SP = Tipo ( 1 d0 ) ! establecer E/S precisión real Privado  Público luban , lusolve Contiene  Subrutina luban ( a , tol , g , h , ip , condinv , detnth ) ! Por Banachiewicz (1938, en adelante B38) El método de descomposición LU calcula tales ! triángulos L=G^T y U=H que el cuadrado B=A^T=G^TH=LU. El pivoteo parcial ! por permutación de columna IP(:) es una adición moderna. ! Dentro del código a, g corresponden a B38 A^T y G^T, de modo que a=gh se cumple. ! ! El uso normal es para el cuadrado A, sin embargo para RHS l ya conocido ! la entrada de (A|l)^T produce (L|y^T)^T donde x en L^Tx=y es la solución de Ax=l. Real ( SP ), Intent ( In ) :: a (:, :) ! matriz de entrada A(m,n), n<=m Real ( SP ), Intent ( In ) :: tol ! tolerancia para pivote cercano a cero Real ( SP ), Intent ( Out ) :: g ( size ( a , dim = 1 ), size ( a , dim = 2 )) ! L(m,n) Real ( SP ), Intent ( Out ) :: h ( size ( a , dim = 2 ), size ( a , dim = 2 )) ! U(n,n) ! nota U columnas están permutadas Real ( SP ), Intent ( Out ) :: condinv ! 1/cond(A), 0 para A singular Real ( SP ), Intent ( Out ) :: detnth ! signo*Abs(det(A))**(1/n) Entero , Intent ( Out ) :: ip( tamaño ( a , dim = 2 )) ! permutación de columnas ! Entero :: k , n , j , l , isig Real ( SP ) :: tol0 , pivmax , pivmin , piv ! n = tamaño ( a , dim = 2 ) tol0 = Máx ( tol , 3._SP * épsilon ( tol0 )) ! usar predeterminado para tol=0 ! ! Se permiten rectangulares A y G bajo la condición: Si ( n > tamaño ( a , dim = 1 ) . O . n < 1 ) Detener 91 Para todo ( k = 1 : n ) ip ( k ) = k h = 0._SP g = 0._SP isig = 1 detnth = 0._SP pivmax = Maxval ( Abs ( a ( 1 , :))) pivmin = pivmax ! Hacer k = 1 , n ! Banachiewicz (1938) Eq. (2.3) h ( k , ip ( k :)) = a ( k , ip ( k :)) - Matmul ( g ( k , : k - 1 ), h (: k - 1 , ip ( k :))) ! ! Encontrar el pivote de fila j = ( Maxloc ( Abs ( h ( k , ip ( k :))), dim = 1 ) + k - 1) Si ( j /= k ) Entonces ! Intercambiar columnas j y k isig = - isig ! Cambiar signo de Det(A) debido a la permutación l = ip ( k ) ip ( k ) = ip ( j ) ip ( j ) = l Fin Si piv = Abs ( h ( k , ip ( k ))) pivmax = Max ( piv , pivmax ) ! Ajustar condinv pivmin = Min ( piv , pivmin ) Si ( piv < tol0 ) Entonces ! matriz singular isig = 0 pivmax = 1._SP Salir  Si no ! Contabilizar la contribución del pivote al signo y valor de Det(A) Si ( h ( k , ip ( k )) < 0._SP ) isig = - isig detnth = detnth + Log ( piv ) Fin Si ! ! Ecuación de Banachiewicz (1938) transpuesta. (2.4) g ( k + 1 :, k ) = ( a ( k + 1 :, ip ( k )) - & Matmul ( g ( k + 1 :, : k - 1 ), h (: k - 1 , ip ( k )))) / h ( k , ip ( k )) g ( k , k ) = 1._SP Fin del bucle ! detnth = isig * Exp ( detnth / n ) condinv = Abs (isig ) * pivmin / pivmax ! Prueba para cuadrado A(n,n) descomentando abajo ! Imprimir *, '|AQ-LU| ',Maxval (Abs(a(:,ip(:))-Matmul(g, h(:,ip(:))))) Fin Subrutina luban Subrutina lusolve ( l , u , ip , x ) ! Resuelve el sistema Ax=b usando factores triangulares LU=A Real ( SP ), Intent ( In ) :: l (:, :) ! matriz triangular inferior L(n,n) Real ( SP ), Intent ( In ) :: u (:, :) ! matriz triangular superior U(n,n) Integer , Intent ( In ) :: ip (:) ! permutación de columnas IP(n) Real ( SP ), Intent ( InOut ) :: x (:, :) ! Entrada: m conjuntos de RHS B(n,m), ! Salida: los conjuntos correspondientes de incógnitas X(n,m) Entero :: n , m , i , j n = tamaño ( ip ) m = tamaño ( x , dim = 2 ) Si ( n < 1. O . m < 1. O . Cualquier ([ n , n ] /= forma ( l )). O . Cualquier ( forma ( l ) /= forma ( u )). O . & n /= tamaño ( x , dim = 1 )) Detener 91 Hacer i = 1 , m Hacer j = 1 , n x ( j , i ) = x ( j , i ) - producto_punto ( x(: j - 1 , i ), l ( j ,: j - 1 )) Fin del bucle  Hacer j = n , 1 , - 1 x ( j , i ) = ( x ( j , i ) - producto_punto ( x ( j + 1 :, i ), u ( j , ip ( j + 1 :)))) / & u ( j , ip ( j )) Fin del  bucle Fin del bucle Fin de  la subrutina lusolve Fin del módulo mlu

Ejemplo de código C

/* ENTRADA: A - matriz de punteros a filas de una matriz cuadrada de dimensión N * Tol - pequeño número de tolerancia para detectar fallos cuando la matriz está cerca de la degeneración * 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 ; double maxA , * ptr , absA ;para ( i = 0 ; i <= N ; i ++ ) P [ i ] = i ; //Matriz de permutación unitaria, P[N] inicializada con Npara ( 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 ; // fallo, la matriz es degeneradaif ( 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 el 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 ]; } }return 1 ; //descomposición realizada }/* ENTRADA: A, P rellenas en LUPDecompose; b - vector del lado 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 rellenas en 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 ;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 rellenas en LUPDecompose; N - 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 ];para ( int i = 1 ; i < N ; i ++ ) det *= A [ i ][ i ];devolver ( P [ N ] - N ) % 2 == 0 ? det : - det ; }

Ejemplo de código C#

public class SystemOfLinearEquations { public double [] SolveUsingLU ( double [,] matrix , double [] rightPart , int n ) { // descomposición de la matriz double [,] lu = new double [ n , n ]; double sum = 0 ; for ( int i = 0 ; i < n ; i ++ ) { for ( int j = i ; j < n ; j ++ ) { sum = 0 ; for ( int k = 0 ; k < i ; k ++ ) sum += lu [ i , k ] * lu [ k , j ]; lu [ i , j ] = matrix [ i , j ] - sum ; } for ( int j = i + 1 ; j < n ; j ++ ) { sum = 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 // encontrar solución de Ly = b double [] y = new double [ n ]; for ( int i = 0 ; i < n ; i ++ ) { sum = 0 ; for ( int k = 0 ; k < i ; k ++ ) sum += lu [ i , k ] * y [ k ]; y [ i ] = rightPart [ i ] - sum ; } // encontrar solución de Ux = y double [] x = new double [ n ]; for ( int i = n - 1 ; i >= 0 ; i -- ) { sum = 0 ; for ( int k = i + 1 ; k < n ; k ++ ) sum += lu [ i , k ] * x [ k ]; x [ i ] = ( 1 / lu [ i , i ]) * ( y [ i ] - sum ); } return x ; } }

Ejemplo de código MATLAB

función LU = LUDecompDoolittle ( A ) n = longitud ( A ); LU = A ; para k = 2 : n para i = 1 : k - 1 lamda = LU ( k , i ) / LU ( i , i ); LU ( k , i ) = lamda ; LU ( k , i + 1 : n ) = LU ( k , i + 1 : n ) - LU ( i , i + 1 : n ) * lamda ; fin fin finfunción x = SolveLinearSystem ( LU, B ) n = length ( LU ); y = zeros ( size ( B )); % encontrar solución de Ly = B para i = 1 : n y ( i ,:) = B ( i ,:) - LU ( i , 1 : i ) * y ( 1 : i ,:); fin % encontrar solución de Ux = y x = zeros ( size ( B )); para i = n :( - 1 ): 1 x ( i ,:) = ( y ( i ,:) - LU ( i ,( i + 1 ): n ) * x (( i + 1 ): n ,:)) / LU ( i , i ); fin finA = [ 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

Notas

↑ Schwarzenberg-Czerny, A. (1995). "Sobre la factorización de matrices y la solución eficiente de mínimos cuadrados" . Astronomy and Astrophysics Supplement Series . 110 : 405. Bibcode : 1995A & AS..110..405S .
↑ Cormen et al. (2009) , pág. 819 , 28.1: Resolución de sistemas de ecuaciones lineales.
1 2 Okunev y Johnson (1997) , Corolario 3 .
↑ Trefethen y Bau (1997) , pág. 166.
↑ Trefethen y Bau (1997) , pág. 161.
↑ Banachiewicz (1938) ; Lay, Lay y McDonald (2021) , pág. 133 , 2.5: Factorizaciones matriciales.
↑ Rigotti (2001) , Principal Menor Principal.
1 2 Horn y Johnson (1985) , Corolario 3.5.5
↑ Horn y Johnson (1985) , Teorema 3.5.2.
↑ Nhiayi, Ly; Phan-Yamada, Tuyetdong (2021). "Examinando la posible descomposición LU". North American GeoGebra Journal . 9 (1).
↑ Okunev y Johnson (1997) .
↑ Householder (1975) .
↑ Golub & Van Loan (1996) , págs.112 , 119.
^ 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 .
↑ Bunch y Hopcroft (1974) .
↑ Trefethen y Bau (1997) , pág. 152.
↑ Préstamo Golub y Van (1996) , pág. 121.
↑ Ralston (1965) .
↑ Hart (2011) .
1 2 Banachiewicz (1938) .
↑ Dwyer (1951) .

Referencias

Banachiewicz, T. (1938), "Méthode de résolution numérique des équations linéaires ..." (PDF) , Bull. Interno. De l'Acad. Polonesa, Serie A. Sc. Matemáticas. : 393– 404.
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), Matrix Computations (3.ª ed.), Baltimore: Johns Hopkins, ISBN 978-0-8018-5414-9.
Hart, Roger (2011), Las raíces chinas del álgebra lineal , Baltimore: Johns Hopkins, ISBN 978-0801897559.
Horn, Roger A.; Johnson, Charles R. (1985), Análisis matricial , Cambridge University Press, ISBN 978-0-521-38632-6Véase la Sección 3.5. N − 1
Householder, Alston S. (1975), La teoría de las 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 ed.), 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.
Ralston, Anthony (1965), Un primer curso de análisis numérico , Nueva York: McGraw-Hill, Inc., ISBN 978-0-070-51157-6.
Press, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007), "Sección 2.3" , Numerical Recipes: The Art of Scientific Computing (3.ª ed.), Nueva York: Cambridge University Press, ISBN 978-0-521-88068-8.
Trefethen, Lloyd N.; Bau, David (1997), Álgebra lineal numérica , Filadelfia: Society for Industrial and Applied Mathematics, ISBN 978-0-89871-361-9.
Rigotti, Luca (2001), ECON 2001 - Introducción a los métodos matemáticos, Lección 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 de MATLAB.

Código informático

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

Recursos en línea

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

[Schwarzenberg-1] Schwarzenberg-Czerny, A. (1995). "Sobre la factorización de matrices y la solución eficiente de mínimos cuadrados" . Astronomy and Astrophysics Supplement Series . 110 : 405. Bibcode : 1995A & AS..110..405S .

[FOOTNOTECormenLeisersonRivestStein2009[httpsarchiveorgdetailsintroduction-to-algorithms-third-edition-2009page819_819]28.1:_Solving_systems_of_linear_equations-2] Cormen et al. (2009) , pág. 819 , 28.1: Resolución de sistemas de ecuaciones lineales.

[okunev-cor3-3] 1 2 Okunev y Johnson (1997) , Corolario 3 .

[FOOTNOTETrefethenBau1997166-4] Trefethen y Bau (1997) , pág. 166.

[FOOTNOTETrefethenBau1997161-5] Trefethen y Bau (1997) , pág. 161.

[FOOTNOTEBanachiewicz1938LayLayMcDonald2021[httpsarchiveorgdetailslinear-algebra-its-applications-6thpage133_133]2.5:_Matrix_Factorizations-6] Banachiewicz (1938) ; Lay, Lay y McDonald (2021) , pág. 133 , 2.5: Factorizaciones matriciales.

[FOOTNOTERigotti2001Leading_Principal_Minor-7] Rigotti (2001) , Principal Menor Principal.

[horn-cor355-8] 1 2 Horn y Johnson (1985) , Corolario 3.5.5

[FOOTNOTEHornJohnson1985Theorem_3.5.2-9] Horn y Johnson (1985) , Teorema 3.5.2.

[10] Nhiayi, Ly; Phan-Yamada, Tuyetdong (2021). "Examinando la posible descomposición LU". North American GeoGebra Journal . 9 (1).

[FOOTNOTEOkunevJohnson1997-11] Okunev y Johnson (1997) .

[FOOTNOTEHouseholder1975-12] Householder (1975) .

[FOOTNOTEGolubVan_Loan1996112,_119-13] Golub & Van Loan (1996) , págs.112 , 119.

[14] 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 .

[FOOTNOTEBunchHopcroft1974-15] Bunch y Hopcroft (1974) .

[FOOTNOTETrefethenBau1997152-16] Trefethen y Bau (1997) , pág. 152.

[FOOTNOTEGolubVan_Loan1996121-17] Préstamo Golub y Van (1996) , pág. 121.

[FOOTNOTERalston1965-18] Ralston (1965) .

[FOOTNOTEHart2011-19] Hart (2011) .

[FOOTNOTEBanachiewicz1938-20] 1 2 Banachiewicz (1938) .

[FOOTNOTEDwyer1951-21] Dwyer (1951) .

[

[ 3 ]

[ 4 ]

de

[ 7 ]

[ 8 ]

[ 9 ]

[ 10 ]

[ 11 ]

[ 12 ]

el

[ 14 ]

[ 15 ]

[

X.

[ 18 ]

[ 19 ]

[ 20 ]

[ 21 ]

Definiciones

LU mediante recursión

Factorización LU con pivoteo parcial

Factorización LU con pivoteo completo

Descomposición diagonal inferior-superior (LDU)

Matrices rectangulares

Ejemplo

Existencia y singularidad

Matrices cuadradas

Matrices simétricas definidas positivas

Matrices generales

Algoritmos

Fórmula cerrada

Utilizando la eliminación gaussiana

Explicación generalizada

Notación

Procedimiento

Ejemplo

Relaciones cuando no se intercambian filas

Descomposición de LU Banachiewicz

Descomposición de LU Crout

Algoritmo aleatorio

Complejidad teórica

Descomposición de matrices dispersas

Aplicaciones

Resolución de ecuaciones lineales

Invertir una matriz

Calcular el determinante

Historia

Ejemplos de código

Ejemplo de código Fortran90

Ejemplo de código C

Ejemplo de código C#

Ejemplo de código MATLAB

Véase también

Notas

Referencias

Enlaces externos

Articulos relacionados