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
Si $A \in \mathbb{R}^{m \times n}$ con $m >= n$ entonces existen $Q \in \mathbb{R}^{m \times m}$ ortogonal y $R \in \mathbb{R}^{m \times n}$ triangular superior tales que $A=QR$:
Comentario: la factorización $QR$ puede escribirse en una forma "delgada" obteniéndose la factorización thin QR:
En este caso $Q_1$ no se llama ortogonal sino con columnas ortonormales, $R_1$ es triangular superior. Si el rank(A) es igual a $n$ (matriz $A$ se nombra de rank completo o full rank) y $R_1$ tiene entradas en la diagonal positivas, la factorización thin QR es única.
Un conjunto de vectores $\{x_1, \dots, x_p\}$ en $\mathbb{R}^m$ ($x_i \in \mathbb{R}^m$)es ortogonal si $x_i^Tx_j=0$ $\forall i\neq j$. Por ejemplo, para un conjunto de $2$ vectores $x_1,x_2$ en $\mathbb{R}^3$ esto se visualiza:
Comentarios:
Si el conjunto $\{x_1,\dots,x_p\}$ en $\mathbb{R}^m$ satisface $x_i^Tx_j= \delta_{ij}= \begin{cases} 1 &\text{ si } i=j,\\ 0 &\text{ si } i\neq j \end{cases}$ se le nombra conjunto ortonormal, esto es, constituye un conjunto ortogonal y cada elemento del conjunto tiene norma Euclidiana igual a $1$: $||x_i||_2 = 1, \forall i=1,\dots,p$.
Si definimos a la matriz $X$ con columnas dadas por cada uno de los vectores del conjunto $\{x_1,\dots, x_p\}$: $X=(x_1, \dots , x_p) \in \mathbb{R}^{m \times p}$ entonces la propiedad de que cada par de columnas satisfaga $x_i^Tx_j=\delta_{ij}$ se puede escribir en notación matricial como $X^TX = I_p$ con $I_p$ la matriz identidad de tamaño $p$ o bien $XX^T=I_m$. A la matriz $X$ se le nombra matriz con columnas ortonormales.
Si cada $x_i$ está en $\mathbb{R}^p$ (en lugar de $\mathbb{R}^m$) entonces construímos a la matriz $X$ como el punto anterior con la diferencia que $X \in \mathbb{R}^{p \times p}$. En este caso $X$ se le nombra matriz ortogonal.
La factorización $QR$ es posible calcularla con las reflexiones de Householder las cuales son matrices simétricas, ortogonales, $R^TR = R^2 = I_m$, $R^{-1}=R$, $det(R)=-1$. Una reflexión de Householder se puede construir a partir de un vector $v \neq 0$ definiendo:
$$R = I_m-\beta v v^T$$con $v \in \mathbb{R}^m - \{0\}$ y $\beta = \frac{2}{v^Tv}$. El vector $v$ se llama vector de Householder. La multiplicación $Rx$ representa la reflexión del vector $x \in \mathbb{R}^m$ a través del hiperplano $span\{v\}^\perp$*
*$span\{v\}$ es el conjunto generado por $v$. Se define como el conjunto de combinaciones lineales de $v$: $span\{v\} = \left \{\displaystyle \sum_{i=1}^m k_i v_i | k_i \in \mathbb{R} \forall i =1,\dots,m \right \}$.
Un dibujo que representa lo anterior es el siguiente en el que se utiliza $u \in \mathbb{R}^m - \{0\}$ , $||u||_2 = 1$ y $R=I_m-2 u u^T$, el reflector elemental alrededor de $u^\perp$:
Comentario: en el dibujo anterior se utiliza el proyector ortogonal elemental sobre el complemento ortogonal de u: $u^\perp = \{x \in \mathbb{R}^m| u^Tx=0\}$ definido como: $P=I_m- u u^T$ y $Px$ es la proyección ortogonal de $x$ sobre $u^\perp$ . Los proyectores ortogonales elementales no son matrices ortogonales, son singulares, son simétricas y $P^2=P$. El proyector ortogonal elemental de $x$ sobre $u^\perp$ tienen $rank$ igual a $m-1$ y el proyector ortogonal de $x$ sobre $span\{u\}$ definido por $I_m-P=uu^T$ tienen $rank$ igual a $1$.
Las reflexiones de Householder pueden utilizarse para hacer ceros por debajo de una entrada de un vector. Por ejemplo:
In [1]:
import numpy as np
In [2]:
x = np.array([1,2,3])
In [3]:
x
Out[3]:
Utilizamos la definición $v=x-||x||_2e_1$ con $e_1=(1,0,0)^T$ vector canónico para construir al vector de Householder:
In [4]:
v = x-np.linalg.norm(x)*np.array([1,0,0]) #uno en la primer entrada pues se desea
#hacer ceros en las entradas restantes
In [5]:
beta = 2/v.dot(v)
Hacemos ceros por debajo de la primera entrada de $x$ haciendo la multiplicación matriz-vector $Rx$. Pero para aplicar $Rx$ no construímos $R = I_3 -\beta vv^T$, en lugar de eso hacemos $Rx = x - \beta v v^Tx$:
In [6]:
x-beta*v*(v.dot(x))
Out[6]:
obsérvese que el resultado de $Rx$ es $(||x||_2,0,0)^T$:
In [14]:
np.linalg.norm(x)
Out[14]:
Sólo para mostrar que $Rx$ construyendo $R$ es equivalente a lo anterior, se construye a continuación $R$:
In [17]:
R = np.eye(3)-beta*np.outer(v,np.transpose(v))
In [18]:
R
Out[18]:
In [19]:
R@x
Out[19]:
Comentario:
In [22]:
v = x+np.linalg.norm(x)*np.array([1,0,0]) #uno en la primer entrada pues se desea
#hacer ceros en las entradas restantes
In [24]:
beta = 2/v.dot(v)
In [25]:
x-beta*v*(v.dot(x))
Out[25]:
¿Cuál definición del vector de Householder usar?
La respuesta a la pregunta tiene que ver con que en cualquiera de las dos definiciones del vector de Householder $v=x \pm ||x||_2 e_1$, la multiplicación $Rx$ refleja $x$ en el primer eje coordenado:
El vector $v^+ = - u_0^+ = x-||x||_2e_1$ refleja $x$ respecto al subespacio $H^+$. El vector $v^- = -u_0^- = x+||x||_2e_1$ refleja $x$ respecto al subespacio $H^-$. Para disminuir los errores por redondeo y evitar el problema de cancelación en la aritmética de punto flotante (ver 1.2.Sistema_de_punto_flotante) se utiliza $v = x+signo(x_1)||x||_2e_1$ donde $signo(x_1) = \begin{cases} 1 &\text{ si } x_1 \geq 0 ,\\ -1 &\text{ si } x_1 < 0 \end{cases}$. La idea de la definción anterior con la función $signo(\cdot)$ es que la reflexión (en el dibujo anterior $-||x||_2e_1$ o $||x||_2e_1$) sea lo más alejada posible de $x$. En el dibujo anterior como $x_1, x_2>0$ entonces se refleja respecto al subespacio $H^-$ quedando su reflexión igual a $-||x||_2e_1$.
Comentarios:
Algoritmo para el cálculo del vector de Householder, nombre de función: $house(\cdot)$
Entrada: $x \in \mathbb{R}^m$.
Salida: $v \in \mathbb{R}^m$ con $v_1 =1$ y $\beta \in \mathbb{R}$ tal que $R=I_m-\beta vv^T$ es ortogonal con $Rx = ||x||_2e_1$.
Definir $m$ como la longitud de $x$.
Definir $norm\text{_2_m} = x(2:m)^Tx(2:m)$.
Definir $v = \left[\begin{array}{c} 1\\ x(2:m) \end{array} \right]$
//Checamos con las siguientes condiciones si $x$ es un múltiplo del vector canónico $e_1$ y el $signo(x_1)$
Si $norm\text{_2_m} = 0$, $x_1 \geq 0$ hacer $\beta=0$. En otro caso:
Si $norm\text{_2_m} = 0$, $x_1 <0$ hacer $\beta=2$. #la referencia Golub, Van Loan tiene un error... es beta=2 no beta=-2, corroborar con 3a edición del libro.
En otro caso:
$norm\text{_x} = \sqrt{x_1^2 + norm\text{_2_m}}$
Si $x_1 \leq 0$ hacer $v_1 = x_1-norm\text{_x}$. En otro caso: $v_1 = - \frac{norm\text{_2_m}}{x_1+norm\text{_x}}$
Definir $\beta = \frac{2v_1^2}{norm\text{_2_m}+v_1^2}$
Definir $v = \frac{v}{v_1}$
Algoritmo para calcular la factorización $QR$ con reflexiones de Householder
Entrada: $A \in \mathbb{R}^{m \times n}$ con $m \geq n$.
Salida: matriz $A$ sobreescrita con la parte esencial del vector de Householder y las entradas de la matriz $R$. Ver comentarios de abajo.
Para j de 1 a n
Calcular $v,\beta$ llamando a la función $house(A(j:m,j))$
$A(j:m, j:n) = A(j:m,j:n)-(\beta v)(v^TA(j:m,j:n))$
Si $j < m$ definir $A(j+1:m,j) = v(2:m-j+1)$.
Comentarios:
Si $v^{(j)} = [0, \dots ,0,1,v_{j+1}^{(j)}, \dots,v_m^{(j)}]^T$ es el $j$-ésimo vector de Householder para $j=1,\dots,m$, y por ejemplo $m=5$, $n=4$ entonces con el algoritmo anterior se obtiene:
para $j=1,\dots,m$.
El almacenamiento de los vectores $v^{(1)}, v^{(2)}, \dots, v^{(n)}$ y las $\beta_j$'s correspondientes se le nombra factor form representation de $Q$.
Típicamente no es necesario construír a la matriz $Q$ en la factorización $QR$. Por ejemplo para realizar la multiplicación: $Q^TC$ con $C \in \mathbb{R}^{m \times p}$ se realiza con:
Para $j$ de $1$ a $n$
$v(j:m) = \left[\begin{array}{c} 1\\ A(j+1:m,j) \end{array} \right]$
*$\beta_j = \frac{2}{1+||A(j+1:m,j)||_2^2}$
*$C(j:m,:) = C(j:m,:) - (\beta_jv(j:m))(v(j:m)^TC(j:m,:))$
¿Qué es lo que está realizando el algoritmo para calcular la factorización QR?
Está aplicando las transformaciones de Householder a la matriz $A$ por la izquierda. Si $Q_j = I_m - \beta_j v^{(j)}v^{T(j)}$ es la $j$-ésima transformación de Householder y $v^{(j)}$ es el $j$-ésimo vector de Householder, entonces en el ciclo for del algoritmo se realiza el producto $Q^TA$ con $Q=Q_1Q_2 \cdots Q_n$. En cada paso del ciclo for se hacen ceros por debajo de una entrada de la matriz:
En el dibujo anterior se muestran los ceros en la matriz $A$ después de aplicar $Q_1, Q_2$ a la matriz, lo cual se representa con el producto: $Q_2Q_1A$. La matriz $Q_3$ se aplicará a continuación y modificará las entradas en rojo. Finalmente $Q^TA = R$ y como $Q$ es ortogonal se tendrá la igualdad: $A=QR$.
Paso 1: encontrar factores $Q,R$ tales que $A=QR$. La $Q$ se puede obtener con la factor form representation, no es necesario calcularla pues para el paso $2$ se puede hacer uso del algoritmo de multiplicación $Q^Tb$ descrito en esta nota.
Paso 2: resolver con el método de sustitución hacia atrás el sistema triangular superior $Rx=Q^Tb$.
Referencias:
G. H. Golub, C. F. Van Loan, Matrix Computations, John Hopkins University Press, 2013.
Ver 0_definiciones_generales para algunas definiciones del álgebra lineal como ortogonalidad, matriz de reflexión de Householder y matrices de rotación.