34.3: Enfoque en LU

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:

1. Resolver\(Ly=b\)\(y\) por sustitución hacia adelante; y luego
2. 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\).

## Import all necessary packages
%matplotlib inline
import scipy.sparse as sparse #this helps to speed up the algorithms, but you will not use it. 
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym
sym.init_printing(use_unicode=True)

## These will allow us to see a cool simulation of the Heat Equation problem (if we compute our answers correctly!)
from matplotlib import animation, rc
from IPython.display import HTML

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\]

Hacer esto

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.

A = np.matrix([[2,1,1,0],[4,3,3,1],[8,7,9,5],[6,7,9,8]]) # Here is A for your convenience
As = sym.Matrix(A)
As

## Type your answer here ##
L1 = np.matrix([[1,,,],[,1,,],[,,1,],[,,,1]])

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.

Hacer esto

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.

## Type your answer here ##

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\]

Hacer esto

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.

## Type your answer here ##
L2 = np.matrix([[1,,,],[,1,,],[,,1,],[,,,1]]) # for your convenience

Hacer esto

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.

## Type your answer here ##

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!

Hacer esto

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.

## Type your answer here ##
L3 = np.matrix([[1,,,],[,1,,],[,,1,],[,,,1]]) # for your convenience

Hacer esto

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!

## Type your answer here ##

¡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\]

Hacer esto

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\)

## Type your answer here ##

Hacer esto

Mirar todas las matrices\(L_i\) y ver las conexiones entre la final\(L\).

print(L1)
print(L2)
print(L3)
print(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.

Hacer esto

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.

## Type your answer here ##

Algoritmo de descomposición general

Hacer esto

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).

## Type your answer here ##
C = np.matrix([[2,1,1,0],[4,3,3,1],[8,7,9,5],[6,7,9,8]]) # to test

def LU_decomp(B):
    n = len(B)
    U = B.copy()
    L = np.identity(n)
    for k in np.arange(0,n-1):
        for j in np.arange(k+1,n):
            L[j,k] =
            U[j,k:n] = U[,:] - L[,]*U[,:]
    return np.mat(L), np.mat(U)

L1,U1 = LU_dec(C) # syntax for returning matrices
np.linalg.norm(L1*U1 - A) # Test: should return 0

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\).

Hacer esto

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.)

## Type your answer here ##
def LU_Axb_solve(A,b):
    L,U = LU_decomp(A)
    n = len(A)
    # Forward-Substitution: Ly = b for y
    y = np.zeros((,))
    for i in np.arange(0,n):
        y[i] = b[i].copy()
        for j in np.arange(0,i):
            y[] = y[] - L[,]*y[]
       
    # Backward-Substitution: Ux = y for x
    x = np.zeros((n,1))
    for i in np.arange(n-1,-1,-1):
        x[] = y[].copy()
        for j in np.arange(n-1,i,-1):
            x[] = x[] - U[,]*x[]
        x[] = x[]/U[,]
    
    return np.mat(x)