34.3: Enfoque en LU
- Page ID
- 115604
\( \newcommand{\vecs}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} } \)
\( \newcommand{\vecd}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash {#1}}} \)
\( \newcommand{\id}{\mathrm{id}}\) \( \newcommand{\Span}{\mathrm{span}}\)
( \newcommand{\kernel}{\mathrm{null}\,}\) \( \newcommand{\range}{\mathrm{range}\,}\)
\( \newcommand{\RealPart}{\mathrm{Re}}\) \( \newcommand{\ImaginaryPart}{\mathrm{Im}}\)
\( \newcommand{\Argument}{\mathrm{Arg}}\) \( \newcommand{\norm}[1]{\| #1 \|}\)
\( \newcommand{\inner}[2]{\langle #1, #2 \rangle}\)
\( \newcommand{\Span}{\mathrm{span}}\)
\( \newcommand{\id}{\mathrm{id}}\)
\( \newcommand{\Span}{\mathrm{span}}\)
\( \newcommand{\kernel}{\mathrm{null}\,}\)
\( \newcommand{\range}{\mathrm{range}\,}\)
\( \newcommand{\RealPart}{\mathrm{Re}}\)
\( \newcommand{\ImaginaryPart}{\mathrm{Im}}\)
\( \newcommand{\Argument}{\mathrm{Arg}}\)
\( \newcommand{\norm}[1]{\| #1 \|}\)
\( \newcommand{\inner}[2]{\langle #1, #2 \rangle}\)
\( \newcommand{\Span}{\mathrm{span}}\) \( \newcommand{\AA}{\unicode[.8,0]{x212B}}\)
\( \newcommand{\vectorA}[1]{\vec{#1}} % arrow\)
\( \newcommand{\vectorAt}[1]{\vec{\text{#1}}} % arrow\)
\( \newcommand{\vectorB}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} } \)
\( \newcommand{\vectorC}[1]{\textbf{#1}} \)
\( \newcommand{\vectorD}[1]{\overrightarrow{#1}} \)
\( \newcommand{\vectorDt}[1]{\overrightarrow{\text{#1}}} \)
\( \newcommand{\vectE}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash{\mathbf {#1}}}} \)
\( \newcommand{\vecs}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} } \)
\( \newcommand{\vecd}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash {#1}}} \)
\(\newcommand{\avec}{\mathbf a}\) \(\newcommand{\bvec}{\mathbf b}\) \(\newcommand{\cvec}{\mathbf c}\) \(\newcommand{\dvec}{\mathbf d}\) \(\newcommand{\dtil}{\widetilde{\mathbf d}}\) \(\newcommand{\evec}{\mathbf e}\) \(\newcommand{\fvec}{\mathbf f}\) \(\newcommand{\nvec}{\mathbf n}\) \(\newcommand{\pvec}{\mathbf p}\) \(\newcommand{\qvec}{\mathbf q}\) \(\newcommand{\svec}{\mathbf s}\) \(\newcommand{\tvec}{\mathbf t}\) \(\newcommand{\uvec}{\mathbf u}\) \(\newcommand{\vvec}{\mathbf v}\) \(\newcommand{\wvec}{\mathbf w}\) \(\newcommand{\xvec}{\mathbf x}\) \(\newcommand{\yvec}{\mathbf y}\) \(\newcommand{\zvec}{\mathbf z}\) \(\newcommand{\rvec}{\mathbf r}\) \(\newcommand{\mvec}{\mathbf m}\) \(\newcommand{\zerovec}{\mathbf 0}\) \(\newcommand{\onevec}{\mathbf 1}\) \(\newcommand{\real}{\mathbb R}\) \(\newcommand{\twovec}[2]{\left[\begin{array}{r}#1 \\ #2 \end{array}\right]}\) \(\newcommand{\ctwovec}[2]{\left[\begin{array}{c}#1 \\ #2 \end{array}\right]}\) \(\newcommand{\threevec}[3]{\left[\begin{array}{r}#1 \\ #2 \\ #3 \end{array}\right]}\) \(\newcommand{\cthreevec}[3]{\left[\begin{array}{c}#1 \\ #2 \\ #3 \end{array}\right]}\) \(\newcommand{\fourvec}[4]{\left[\begin{array}{r}#1 \\ #2 \\ #3 \\ #4 \end{array}\right]}\) \(\newcommand{\cfourvec}[4]{\left[\begin{array}{c}#1 \\ #2 \\ #3 \\ #4 \end{array}\right]}\) \(\newcommand{\fivevec}[5]{\left[\begin{array}{r}#1 \\ #2 \\ #3 \\ #4 \\ #5 \\ \end{array}\right]}\) \(\newcommand{\cfivevec}[5]{\left[\begin{array}{c}#1 \\ #2 \\ #3 \\ #4 \\ #5 \\ \end{array}\right]}\) \(\newcommand{\mattwo}[4]{\left[\begin{array}{rr}#1 \amp #2 \\ #3 \amp #4 \\ \end{array}\right]}\) \(\newcommand{\laspan}[1]{\text{Span}\{#1\}}\) \(\newcommand{\bcal}{\cal B}\) \(\newcommand{\ccal}{\cal C}\) \(\newcommand{\scal}{\cal S}\) \(\newcommand{\wcal}{\cal W}\) \(\newcommand{\ecal}{\cal E}\) \(\newcommand{\coords}[2]{\left\{#1\right\}_{#2}}\) \(\newcommand{\gray}[1]{\color{gray}{#1}}\) \(\newcommand{\lgray}[1]{\color{lightgray}{#1}}\) \(\newcommand{\rank}{\operatorname{rank}}\) \(\newcommand{\row}{\text{Row}}\) \(\newcommand{\col}{\text{Col}}\) \(\renewcommand{\row}{\text{Row}}\) \(\newcommand{\nul}{\text{Nul}}\) \(\newcommand{\var}{\text{Var}}\) \(\newcommand{\corr}{\text{corr}}\) \(\newcommand{\len}[1]{\left|#1\right|}\) \(\newcommand{\bbar}{\overline{\bvec}}\) \(\newcommand{\bhat}{\widehat{\bvec}}\) \(\newcommand{\bperp}{\bvec^\perp}\) \(\newcommand{\xhat}{\widehat{\xvec}}\) \(\newcommand{\vhat}{\widehat{\vvec}}\) \(\newcommand{\uhat}{\widehat{\uvec}}\) \(\newcommand{\what}{\widehat{\wvec}}\) \(\newcommand{\Sighat}{\widehat{\Sigma}}\) \(\newcommand{\lt}{<}\) \(\newcommand{\gt}{>}\) \(\newcommand{\amp}{&}\) \(\definecolor{fillinmathshade}{gray}{0.9}\)En esta asignación crearemos algoritmos que factorizan matrices invertibles (es decir, matrices cuadradas con columnas linealmente independientes)\(A\), en la siguiente descomposición (nota: los términos descomposición y factorización se usan indistintamente en la literatura ):
- Descomposición de LU:\(A=LU\)
Cada matriz en estas descomposiciones es el producto matricial de matrices elementales. Recordemos que una matriz elemental difiere de la matriz de identidad (la matriz cuadrada con\(1\) s en la diagonal y\(0\) s en otra parte) por una operación de fila elemental.
El uso de estas descomposiciones matriciales en Álgebra Lineal Numérica está motivado por la eficiencia computacional o la reducción de la complejidad computacional (recuerde “notación Big-O” y contando los bucles en su algoritmo de multiplicación matricial) y estabilidad numérica . Resolver a nuestro viejo amigo\(Ax=b\) calculando lo inverso de\(A\)\(A^{−1}\), denotado y tomando el producto matriz-vector\(A^{−1}b=x\) es intensivo en recursos computacionales y numéricamente inestable, en general. Si la\(LU\) descomposición de\(A\) existe, entonces será menos costosa y más estable para:
- Resolver\(Ly=b\)\(y\) por sustitución hacia adelante; y luego
- Resolver\(Ux=y\)\(x\) por sustitución hacia atrás
Una nota final para relacionar esta tarea con el inicio del curso: Los algoritmos aquí presentados son de una clase diferente a la del Algoritmo Jacobi y el Algoritmo de Gauss-Siedel. Estos son algoritmos iterativos. Como ya sabe, esto significa que el procedimiento algorítmico se aplica una, dos veces, y así sucesivamente, hasta que la salida esté dentro de una tolerancia deseada, o hasta que el proceso se haya ejecutado un número determinado de veces (por ejemplo, 100 iteraciones).
Eliminación gaussiana y descomposición de LU
Recordemos que la eliminación gaussiana convierte una matriz arbitraria\(A\),, en su forma de escalón de filas. Para nuestros propósitos, supongamos que\(A\) es una matriz cuadrada y, por lo tanto, una\(n \times n\) matriz. Para simplificar nuestras tareas, impongamos aún más la condición que\(A\) es invertible. Así, las columnas de\(A\) son linealmente independientes. Esto significa que la eliminación gaussiana dará como resultado una matriz triangular superior. Denotemos esta matriz\(U\) para triangulares superiores.
Si hubiera una función,\(f\) que pudiera tomar\(A\) y dar salida\(U\), podríamos pensar en la Eliminación Gaussiana como el siguiente proceso:
\[f(A)=U \nonumber\]
Con esta información, ahora podemos reescribir nuestra ecuación desde arriba como:
\[L^{-1}A = U \nonumber\]
Es posible que hayas notado el superíndice en\(L^{−1}\). Esto solo dice que\(L^{−1}\) es lo inverso de alguna matriz\(L\). Y para cualquier matriz invertible,\(L\), tenemos que los productos de la matriz:
\[L^{-1}L = LL^{-1} = I \nonumber\]
Esto es análogo a (para cada número real\(a \neq 0\)):
\[a^{-1}\times a = a\times a^{-1} = 1 \nonumber\]
Usando las reglas de multiplicación matricial, verifica la fórmula anterior calculando lo siguiente:
\ [\ begin {split}
L_ {1} ^ {-1} L_ {1} =
\ left (
\ begin {array} {*5 {c}}
1 & 0 & 0 & 0 & 0 & 0\
-l_ {21} & 1 & 0 & 0 & 0\\
-l_ {31} & 0 & 0 & 0 & 0\\
-l_ {41} & 0 & 0 & 1 & 0\
-l_ {51} & 0 & 0 & 0 & 1\
\ end {array}
\ derecha)\ left (
\ begin {array} {*5 {c}}
1 & 0 & 0 & 0 & 0\\
l_ {21} & 1 & 0 & 0 & 0\\
l_ {31} & 0 & 1 & 0 & 0\\
l_ {41} & 0 & 0 & 0 & 1 & 0\
l_ {51} & 0 & 0 & 0 & 1\\
\ end {array}
\ right)
=
\ left (
\ begin {array} {*5 {c}}
1 & 0 & 0 & 0 & 0 & 0 &
0 & 1 & 0 & 0 & 0 &
0 & 0\\ 0 & 0 & 0 \\
0 & 0 & 0 & 1 & 0\\
0 & 0 & 0 & 0 & 0 & 1\
\\ end {array}
\ right)
= I
\ end {split}\ nonumber\]
\ [\ begin {split}
L_ {2} ^ {-1} L_ {2} =
\ left (
\ begin {array} {*5 {c}}
1 & 0 & 0 & 0 & 0 & 0 &
0 & 1 & 0 & 0 & 0 & 0\\
0 & -l_ {32} & 1 & 0 & 0\\
0 & -l_ {42} & 0 & 1 & 0\\
0 & -l_ {52} & 0 & 0 & 1\
\ end {array}
\ derecha)\ left (
\ begin {array} {*5 {c}}
1 & 0 & 0 & 0 & 0\ \
0 & 1 & 0 & 0 & 0\\
0 & l_ {32} & 1 & 0 & 0\\
0 & l_ {42} & 0 & 0 & 1 & 0\\
0 & l_ {52} & 0 & 0 & 1\\
\ end {array}
\ derecha)
=
\ izquierda (
\ begin {array} {*5 {c}}
1 & 0 & 0 & 0 & 0 & 0 & 1 &
0 & 0 & 0 & 0 & 0\\
0 & 0 & 0 & 0 & 0 &
0 & 0 & 0 & 1 & 0\\
0 & 0 & 0 & 0 & 0 & 1\
\\ end {array}
\ right)
= I
\ end {split}\ nonumber\]
\ [\ begin {split}
L_ {4} ^ {-1} L_ {4} =
\ left (
\ begin {array} {*5 {c}}
1 & 0 & 0 & 0 & 0 &
0\\ 0 & 1 & 0 & 0 &
0 & 0 & 0 & 0 & 0 & 0\
0 & 0 & 0 & 1 & 0\\
0 & 0 & 0 & -l_ {54} & 1\
\ end {array}
\ derecha)\ left (
\ begin {array} {*5 {c}}
1 & 0 & 0 & 0 & 0 & 0\\
0 & 1 & 0 & 0 & 0\\
0 & 0 & 1 & 0 & 0 &
0\\ 0 & 0 & 0 &
0 & 0 & 0 & l_ {54} & 1\\
\ end {array}
\ derecha)
=
\ izquierda (
\ begin {array} {*5 {c}}
1 & 0 & 0 & 0 & 0 & 0 & 1 &
0 & 0 & 0 & 0 & 0\\
0 & 0 & 0 & 0 & 0 &
0 & 0 & 0 & 1 & 0\\
0 & 0 & 0 & 0 & 0 & 1\
\\ end {array}
\ right)
= I
\ end {split}\ nonumber\]
Para entender\(L^{−1}\) más profundamente, volvamos nuestra atención hacia la eliminación gaussiana por un momento. Tomemos como un hecho que, para una\(n \times n\) matriz “suficientemente agradable”\(A\), la matriz\(L^{−1}\) que toma\(A\) a su forma de escalón superior triangular o fila,\(U\), tiene la estructura:
\[L^{-1} = L_{n-1}L_{n-2} \ldots L_{2}L_{1} \nonumber\]
Cada una de las\(L_i\) está arriba es una matriz elemental que pone a cero las entradas subdiagonales de la\(i^{th}\) columna de\(A\). Este es el\(i^{th}\) paso de Eliminación Gaussiana aplicado a toda la\(i^{th}\) columna de A debajo del elemento\(i^{th}\) diagonal.
Demostremos esto por cómputo de\(L_i\) para una matriz “agradable”\(A\).
Eliminación Gaussiana por Matrices Elementales,\(L_i\)
Dejar\(A\) ser la siguiente matriz:
\ [\ begin {split} A =
\ begin {bmatrix}
2 & 1 & 1 & 0\\
4 & 3 & 3 & 1\\
8 & 7 & 9 & 5\\
6 & 7 & 9 & 8\
\ end {bmatrix}
\ end {split}\ nonumber\]
Crear una matriz triangular inferior\(4 \times 4\) unitaria,\(L_1\) que elimina todas las entradas subdiagonales de la primera columna de\(A\). Muestre la matriz\(L_1\) usando SympY.
Ahora deberíamos tener lo siguiente:
\ [\ begin {split} L_ {1} A = A^ {(1)} =
\ begin {bmatrix}
2 & 1 & 1 & 0\\
0 & u_ {22} & u_ {23} & u_ {24}\\
0 & x & x & x\
0 & x & x & x & x\
\ end {bmatrix}
=
\ begin {bmatrix}
u_ {11} & u_ {12} & u_ {13} & u_ {14}\\
0 & u_ {22} & u_ {23} & u_ {24}\\
0 & x & x & x & x\
0 & x & x & x & x\
\ fin { bmatrix}
\ end {split}\ nonumber\]
Dado que nuestra primera fila se mantuvo sin cambios, sabemos que nuestras\(u_{1i}\) (las entradas de la primera fila de\(U\)) ahora están determinadas. De igual manera, tenemos que también se determinan las\(u_{2i}\) (las entradas de la segunda fila de\(U\)). Los\(x\) elementos son elementos que han cambiado, pero que aún no están en su forma final. Nota: Tu\(u_{ij}\) será entero, o integer (\(\mathbb{Z}\)), números.
Multiplicar\(L_1\) a la izquierda\(A\) por para confirmar que todas las entradas subdiagonales de la primera columna de\(A^{(1)}\) son cero. Mostrar el resultado a través de SympY.
Nuestro siguiente paso será eliminar todas las entradas distintas de cero de la segunda columna de\(A^{(1)}=L_{1}A\) por multiplicación izquierda de\(L_2\). Esto debería producir:
\ [\ begin {split}\ begin {align} A^ {(2)} &= L_ {2} A^ {(1)}\ nonumber\\
&= L_ {2} L_ {1} A\ nonumber\\
&=
\ begin {bmatrix}
u_ {11} & u_ {12} & u_ {13} & u_ {14}\\
0 & u_ _ {22} & u_ {23} & u_ {24}\\
0 & 0 & u_ {33} & u_ {34}\\
0 & 0 & x & x\
\ final {bmatrix}\ nonumber
\ end {align}
\ end {split}\ nonumber\]
Crear una matriz triangular inferior\(4 \times 4\) unitaria,\(L_2\) que elimine todas las entradas subdiagonales de la segunda columna de\(A^{(1)}\) rendimiento\(A^{(2)}\) como arriba. Muestre la matriz\(L_2\) usando SympY.
Multiplicar\(L_2\) a la izquierda\(A^{(1)}\) por para confirmar que todas las entradas subdiagonales de la columna 2 de\(A^{(2)}\) son cero. Mostrar el resultado a través de SympY. Nota: Tu\(u_{ij}\) será entero, o entero (\(\mathbb{Z}\)), números.
Ahora deberíamos tener:
\ [\ begin {split}
\ begin {align} A^ {(2)} &= L_ {2} A^ {(1)}\ nonumber\\
&= L_ {2} L_ {1} A\ nonumber\\
&=
\ begin {bmatrix}
u_ {11} & u_ {12} & u_ {13} & u_ {14}\\
0 & u_ _ {22} & u_ {23} & u_ { 24}\\
0 & 0 & u_ {33} & u_ {34}\\
0 & 0 & x & x\
\ final {bmatrix}\ nonumber
\ end {align}
\ end {split}\ nonumber\]
Ahora queremos construir la matriz final\(L_3\) que llevará nuestra matriz\(A^{(2)}\) a forma triangular superior. Entonces, ¡hagámoslo!
Cree una matriz triangular inferior\(4 \times 4\) unitaria,\(L_3\) que elimine todas las entradas subdiagonales de la tercera columna de\(A^{(2)}\) rendimiento:
\ [
\ begin {align} A^ {(3)} &= L_ {3} A^ {(2)}\ nonumber\\
&= L_ {3} L_ {2} A^ {(1)}\ nonumber\\
&= L_ {3} L_ {2} L_ {1} A\ nonumber\\
&= U\ nonumber
\ end align {}
\ nonumber\]
Muestre la matriz\(L_3\) usando SympY.
Multiplicar\(L_3\) a la izquierda\(A^{(2)}\) por para confirmar que todas las entradas subdiagonales de la columna 3 de\(A^{(3)}\) son cero. Mostrar el resultado a través de SympY. Nota: Tu\(u_{ij}\) será entero, o integer (\(\mathbb{Z}\)), números. Ahora debes notar que\(A^{(3)} = U\) está en forma de escalón de fila, y, por lo tanto,\(U\) es una matriz superior-triangular con\(0\) s por debajo de la diagonal!
¡Felicidades!
Acabas de descomponer tu primera matriz a través del proceso de abajo (y ahora debería tener una matriz,\(U\), que se parece a la de abajo):
\ [\ begin {split}
\ begin {align} L^ {-1} A &= L_ {3} L_ {2} L_ {1} A\ nonumber\\
&= L_ {3} L_ {2} A^ {(1)}\ nonumber\\
&= L_ {3} A^ {(2)}\ nonumber\\
&= A^ {(3)})}\ nonumber\\
&= U\ nonumber\\
&=
\ begin { bmatrix}
2 & 1 & 1 & 0\\
0 & 1 & 1 & 1 & 1\\
0 & 0 & 0 & 2 & 2\\
0 & 0 & 0 & 2
\ end {bmatrix}\ nonumber
\ end {align}
\ end {split}\ nonumber\]
Por último, vamos a generar explícitamente las matrices\(L^{−1}\) y\(L\). Después, muéstralos usando SympY.
Será útil utilizar lo siguiente:
\[\begin{align}L^{-1} &= L_{n-1}L_{n-2} \ldots L_{2}L_{1} \nonumber \end{align} \nonumber\]
y
\ [\ begin {align} L &= (L^ {-1}) ^ {-1}\ nonumber\\
&= (L_ {n-1} L_ {n-2}... L_ {2} L_ {1}) ^ {-1}\ nonumber\\
&= L_ {1} ^ {-1} L_ {2} ^ {-1}... L_ {n-2} ^ -1 {} L_ {n-1} ^ {-1}\ nonumber
\ end {align}
\ nonumber\]
Si estás atascado, consulta el párrafo al inicio de esta sección para obtener la fórmula explícita. Recordar:\(L^{-1}L = LL^{-1} = I\)
Mirar todas las matrices\(L_i\) y ver las conexiones entre la final\(L\).
Para nuestro último poco de diversión de descomposición de LU, confirmemos que tus matrices\(L\) y\(U\) cumplen con la igualdad:
\[A = LU \nonumber\]
En efecto, hay una función en SympY que calculará la descomposición de LU para nosotros.
Ejecute la siguiente función e imprima sus salidas:
\[\texttt{L_actual, U_actual, _ = As.LUdecomposition()} \nonumber\]
Luego, compute:
\[\texttt{L_actual*U_actual - As} \nonumber\]
y confirmar que da salida a la matriz cero.
Algoritmo de descomposición general
Usando el siguiente código de andamiaje, complete el algoritmo de descomposición de LU. (Puede ser útil probar su código en la matriz\(A\) desde arriba).
Resolver\(Ax=b\) vía descomposición LU
Es posible que desee consultar la introducción de esta tarea para obtener una visión general de cómo usar LU Descomposición para resolver\(Ax=b\).
Usando el código de andamiaje a continuación, complete el algoritmo de solucionador LU. El algoritmo debe resolver\(Ly = b\) para\(y\) vía Forward-Substitution y luego\(Ux=y\) para\(x\) por Backward-Substitution. (Puede ser útil probar su código en una matriz\(A\) y vector\(b\) de la tarea 1 u otra fuente.)