Notas para contenedor de docker:
Comando de docker para ejecución de la nota de forma local:
nota: cambiar <ruta a mi directorio> por la ruta de directorio que se desea mapear a /datos dentro del contenedor de docker.
docker run --rm -v <ruta a mi directorio>:/datos --name jupyterlab_numerical -p 8888:8888 -d palmoreck/jupyterlab_numerical:1.1.0password para jupyterlab: qwerty
Detener el contenedor de docker:
docker stop jupyterlab_numericalDocumentación de la imagen de docker palmoreck/jupyterlab_numerical:1.1.0 en liga.
Nota basada en liga
En esta nota suponemos que $A \in \mathbb{R}^{n \times n}$ y $A$ es una matriz con entradas $a_{ij} \in \mathbb{R}^{n \times n} \forall i,j=1,2,\dots,n$.
Esta factorización se puede obtener a partir del método de la eliminación Gaussiana que utiliza transformaciones de Gauss para la fase de eliminación. Ver Gaussian elimination.
La idea de la factorización $LU$ simple es utilizar las transformaciones de Gauss para hacer ceros por debajo de un elemento nombrado pivote. Una transformación de Gauss está definida como $L_k = I_n - \ell_ke_k^T$ con $\ell_k = (0,0,\dots,\ell_{k+1,k},\dots,\ell_{n,k})^T$ y $\ell_{i,k}=\frac{a_{ik}}{a_{kk}} \forall i=k+1,\dots,n$ y $e_k$ es el $k$-ésimo vector canónico: vector con un $1$ en la posición $k$ y ceros en las entradas restantes.
Ejemplo:
Considérese al vector $a=(-2,3,4)^T$. Definir una transformación de Gauss para hacer ceros por debajo de $a_1$ y otra transformación de Gauss para hacer cero la entrada $a_3$
In [1]:
import numpy as np
Para hacer ceros por debajo del pivote $a_1$:
In [2]:
v=np.array([-2,3,4])
In [3]:
pivote = v[0]
In [4]:
l1 = np.array([0,v[1]/pivote, v[2]/pivote])
In [5]:
e1 = np.array([1,0,0])
Obsérvese que por la definición de la transformación de Gauss, no necesitamos construir a la matriz $L_k$
In [6]:
L1v = v-l1*(e1.dot(v))
In [7]:
L1v
Out[7]:
Sólo para mostrar que $L_1 v$ construyendo $L_1$ es equivalente a lo anterior, se construye a continuación $L_1$:
In [8]:
L1 = np.eye(3) - np.outer(l1,e1)
In [9]:
L1
Out[9]:
In [10]:
L1@v
Out[10]:
Para hacer ceros por debajo del pivote $a_2$:
In [11]:
v=np.array([-2,3,4])
In [12]:
pivote = v[1]
In [13]:
l2 = np.array([0,0, v[2]/pivote])
In [14]:
e2 = np.array([0,1,0])
Obsérvese que por la definición de la transformación de Gauss, no necesitamos construir a la matriz $L_k$
In [15]:
L2v = v-l2*(e2.dot(v))
In [16]:
L2v
Out[16]:
Sólo para mostrar que $L_2 v$ construyendo $L_2$ es equivalente a lo anterior, se construye a continuación $L_2$:
In [17]:
L2 = np.eye(3) - np.outer(l2,e2)
In [18]:
L2
Out[18]:
In [19]:
L2@v
Out[19]:
Entre sus ventajas se encuentra el poder usarse para resolver múltiples SEL una vez calculada. Por ejemplo, si queremos resolver $t$ SEL $Ax = b_1, Ax=b_2, \dots, Ax=b_t$ con eliminación Gaussiana resultarían $\mathcal{O}(\frac{2}{3}tn^3)$ flops. Si usamos factorización $LU$ entonces factorizamos $A=LU$ con un costo de $\mathcal{O}(\frac{2}{3}n^3)$ flops y posteriormente resolvemos los $2t$ SEL triangulares con los métodos de sustitución hacia delante y hacia atrás. Lo anterior generarían $\mathcal{O}(\frac{2}{3}n^3) + \mathcal{O}(2tn^2)$ flops que es menos costoso que con eliminación Gaussiana si $t>1$.
A la eliminación Gaussiana y a la factorización LU se les tiene que añadir una búsqueda de elementos con magnitud más grande para tener la propiedad de estabilidad numérica* bajo el redondeo en un SPF (ver 1.2.Sistema_de_punto_flotante, 1.3.Condicion_de_un_problema_y_estabilidad_de_un_algoritmo). Tal búsqueda se puede realizar por renglones, por columnas o en ambas. A la búsqueda por renglones se le llama pivoteo parcial*.
*Aunque al añadir el pivoteo parcial a la eliminación Gaussiana o a la factorización LU, se obtienen métodos que en general son númericamente estables bajo el redondeo, no lo son en la totalidad de ejemplos. Hay ejemplos construídos que destruyen su estabilidad numérica. Sin embargo, tales ejemplos no surgen en la práctica (esto se ha comprobado a través de la historia en el uso de éste método)
*El nombre de pivoteo parcial se utiliza por la búsqueda de pivotes que son los elementos de la matriz por los cuales la eliminación Gaussiana realiza la fase de eliminación. El nombre "parcial" se utiliza para diferenciarla del pivoteo completo en el que se buscan los pivotes tanto en renglones como en columnas. En el pivoteo parcial la búsqueda se hace por renglones.
Teorema
Si $A \in \mathbb{R}^{n \times n}$ es no singular, existen matrices $P, L, U$ de permutación, triangular inferior con $1$'s en la diagonal y triangular superior invertible respectivamente tales que $PA = LU$.
Ejemplo numérico
Sea $$ A= \left[ \begin{array}{ccc} 1 &4 &-2\\ -3 &9 &8\\ 5 &1 &-6\\ \end{array} \right] $$
Calcular la factorización $PA=LU$.
Paso 1
En el primer paso de la factorización $LU$ se busca en la primer columna el pivote que es el elemento con máxima magnitud (si hay empate tomamos el primer índice donde ocurre). En el ejemplo anterior es $|5|=5$ por lo que la primer matriz de permutación $P_1$ se define como:
y se multiplica $P_1A$ resultando en:
La matriz de transformación de Gauss entonces es:
$$ L_1= \left[ \begin{array}{ccc} 1 &0 &0\\ \frac{3}{5} &1 &0\\ -\frac{1}{5} &0 &1\\ \end{array} \right] $$y haciendo la multiplicación $L_1P_1A$ resulta:
Paso 2
La matriz $P_2$ es igual a:
$$ P_2= \left[ \begin{array}{ccc} 1 &0 &0\\ 0 &1 &0\\ 0 &0 &1\\ \end{array} \right] $$La matriz de transformación de Gauss $L_2$ es:
$$ L_2= \left[ \begin{array}{ccc} 1 &0 &0\\ 0 &1 &0\\ 0 &\frac{-19}{48} &1\\ \end{array} \right] $$y el producto $L_2P_2L_1P_1A$ es:
El anterior producto forma a la matriz $U$ de $PA=LU$: $U=L_2P_2L_1P_1A$.
El producto $\hat{L}=P_1L_1^{-1}P_2L_2^{-1}$ ayuda a construír a la $L$ de $PA=LU$:
y la $L$ es: $L=P\hat{L}$:
$$ \left[ \begin{array}{ccc} 1 &0 &0\\ \frac{-3}{5} &1 &0\\ \frac{1}{5} &\frac{19}{48} &1\\ \end{array}\right] $$con $P=P_2P_1$:
Obs: el ejemplo numérico anterior sólo es una guía, no se implementa construyendo las matrices de permutación ni las transformaciones de Gauss. Además obsérvese que se realizan $n-1=2$ pasos.
Algoritmo gaxpy LU con pivoteo parcial
El pivoteo parcial se implementa como una búsqueda por la entrada de máxima magnitud por renglones y un intercambio de renglones con el índice donde ocurre tal máximo.
El algoritmo utiliza un vector $piv \in \mathbb{R}^{n-1}$ (obsérvese que tiene dimensión $n-1$) para codificar a las matrices de permutación que permiten realizar el intercambio de renglones. Una matriz de permutación es aquella que tiene exactamente un $1$ en cada renglón y columna. Por ejemplo:
$$ P= \left[ \begin{array}{cccc} 0 &1 &0 &0\\ 0 &0 &0 &1\\ 0 &0 &1 &0\\ 1 &0 &0 &0 \end{array} \right] $$Si se desea realizar la multiplicación $Px$ con $x=(x_1,x_2,x_3,x_4)^T = (0,-2,3,1)^T$, se codifica con un vector $piv \in \mathbb{R}^{n}$ (obsérvese que tiene dimensión $n$) con entradas iguales a:
$$ piv= \left[ \begin{array}{c} 2\\ 4\\ 3\\ 1 \end{array} \right] $$y la multiplicación $Px$ se realiza con el código:
Para k de 1 a n
intercambiar posición k de x con posición k de pivEn el ejemplo anterior $Px = (x_2,x_4,x_3,x_1)^T$:
$$ Px= \begin{array}{l} \left[ \begin{array}{cc} 0 &1 &0 &0\\ 0 &0 &0 &1\\ 0 &0 &1 &0\\ 1 &0 &0 &0 \end{array} \right] \cdot \left[\begin{array}{c} x_1\\ x_2\\ x_3\\ x_4 \end{array} \right] = \left[\begin{array}{c} x_2\\ x_4\\ x_3\\ x_1 \end{array} \right] = \left[\begin{array}{c} -2\\ 1\\ 3\\ 0 \end{array} \right] \end{array} $$La multiplicación entre matrices de permutación se puede codificar considerando $piv \in \mathbb{R}^{n-1}$ (obsérvese que tiene dimensión $n-1$). Por ejemplo si
$$P_1= \left[\begin{array}{ccc} 0 & 0 & 1\\ 0 & 1 & 0\\ 1 & 0 & 0 \end{array} \right] P_2= \left[\begin{array}{ccc} 1 & 0 & 0\\ 0 & 0 & 1\\ 0 & 1 & 0 \end{array} \right] $$y se desea realizar la multiplicación: $P = P_2P_1$:
$$P=P_2P_1 = \left[\begin{array}{ccc} 1 & 0 & 0\\ 0 & 0 & 1\\ 0 & 1 & 0 \end{array} \right] \left[\begin{array}{ccc} 0 & 0 & 1\\ 0 & 1 & 0\\ 1 & 0 & 0 \end{array} \right] = \left[\begin{array}{ccc} 0 & 0 & 1\\ 1 & 0 & 0\\ 0 & 1 & 0 \end{array} \right] $$entonces $piv=(3,3)^T$ representando la primera posición $3$ un intercambio entre el primer renglón de la matriz identidad con el tercer renglón de la matriz identidad. La segunda posición de $piv$ igual a $3$ representa un intercambio entre el segundo renglón de la matriz permutada anterior con el tercer renglón.
Para observar esto con la notación de matrices considérese a la matriz identidad:
$$ I_3= \left[\begin{array}{ccc} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end{array} \right] $$Primera entrada de $piv$ es $3$ entonces:
$$ I_3= \left[\begin{array}{ccc} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end{array} \right] \underset{3}{\overset{1}{\longleftrightarrow}} P_1= \left[\begin{array}{ccc} 0 & 0 & 1\\ 0 & 1 & 0\\ 1 & 0 & 0 \end{array} \right] $$y segunda entrada de $piv$ es $3$ entonces partiendo de la matriz anterior:
$$ P_1= \left[\begin{array}{ccc} 0 & 0 & 1\\ 0 & 1 & 0\\ 1 & 0 & 0 \end{array} \right] \underset{3}{\overset{2}{\longleftrightarrow}} P= \left[\begin{array}{ccc} 0 & 0 & 1\\ 1 & 0 & 0\\ 0 & 1 & 0\\ \end{array} \right] $$que es el resultado previamente calculado $P_2P_1$.
A continuación se presenta el algoritmo para su cáclulo en una versión gaxpy, ver 3.1.Vectorizacion_BLAS_y_el_uso_del_cache_eficientemente para definición de la operación gaxpy.
Algoritmo (se considera que los índices de los arrays inician en el valor $1$ y no en el $0$)
Entrada: $A \in \mathbb{R}^{n \times n}$.
Salida: factores $L, U \in \mathbb{R}^{n \times n}$, vector $piv \in \mathbb{R}^{n}$. El vector $piv$ se utiliza para construir a la matriz $P$.
Inicializar L con matriz identidad. Inicializar U como una matriz de ceros.
Para $j$ de $1$ a $n$ hacer:
Si $j$ es 1 definir $v = A(:,1)$ # v es la primer columna de $A$.
En otro caso:
Calcular $\hat{a} = P_{j-1} \cdots P_1 A(:,j)$ #aquí se utiliza al vector $piv$ para evitar hacer multiplicaciones de matriz vector, sólo se realizan intercambio de renglones
Resolver $L(1:j-1,1:j-1)z = \hat{a}(1:j-1)$ para $z \in \mathbb{R}^{j-1}$ por sustitución hacia delante
Definir $U(1:j-1,j) = z$
Definir $v(j:n) = \hat{a}(j:n) - L(j:n, 1:j-1)*z$
Determinar $\mu$ con $ j \leq \mu \leq n$ tal que $|v(\mu)| = ||v(j:n)||_\infty$
Definir $piv(j)=\mu$ # piv es un vector de tamaño $n-1$
$v(j) \leftrightarrow v(\mu)$ # $\leftrightarrow$ significa intercambiar $j$-ésima entrada de $v$ por $\mu$-ésima entrada de $v$.
$L(j,1:j-1) \leftrightarrow L(\mu,1:j-1)$ #$\leftrightarrow$ significa intercambio
Definir $U(j,j) = v(j)$.
Si $v(j) \neq 0 $ definir $L(j+1:n,j) = v(j+1:n)/v(j)$
Paso 1: encontrar factores $P,L,U$ tales que $PA=LU$.
Paso 2: resolver con el método de sustitución hacia delante el sistema triangular inferior $Ld=Pb$.
Paso 3: resolver con el método de sustitución hacia atrás el sistema triangular superior $Ux=d$.
Referencias: