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 lo que sigue se supone una matriz cuadrada $A \in \mathbb{R}^{nxn}$.
El número $\lambda$ (real o complejo) se denomina eigenvalor de A si existe $v \in \mathbb{C}^n - \{0\}$ tal que $Av = \lambda v$. El vector $v$ se nombra eigenvector (vector propio o característico) de $A$ correspondiente al eigenvalor $\lambda$.
Obs:
Nota: Si A es simétrica entonces tiene eigenvalores reales y aún más: $A$ tiene eigenvectores reales linealmente independientes y forman un conjunto ortonormal. Entonces $A$ se escribe como un producto de tres matrices nombrado descomposición espectral: $$A = Q \Lambda Q^T$$ donde: $Q$ es una matriz ortogonal cuyas columnas son eigenvectores de $A$ y $\Lambda$ es una matriz diagonal con eigenvalores de $A$.
En lo que sigue se supone $A \in \mathbb{R}^{mxn}$.
El número $\sigma$ se denomina valor singular de $A$ si $\sigma = \sqrt{\lambda_{A^TA}} = \sqrt{\lambda_{AA^T}}$ donde: $\lambda_{A^TA}$ y $\lambda_{AA^T}$ es eigenvalor de $A^TA$ y $AA^T$ respectivamente .
Obs: la definición se realiza sobre $A^TA$ o $AA^T$ pues éstas matrices tienen el mismo espectro y además sus eigenvalores son reales y positivos por lo que $\sigma \in \mathbb{R}$ y de hecho $\sigma \geq 0$ (la raíz cuadrada se calcula para un eigenvalor no negativo).
Asociado con cada valor singular $\sigma$ existen vectores singulares $u,v$ que cumplen con la igualdad: $$Av = \sigma u .$$ Al vector $u$ se le nombra vector singular izquierdo y al vector $v$ se le nombra vector singular derecho.
Si $A \in \mathbb{R}^{mxn}$ entonces existen $U \in \mathbb{R}^{mxm}, V \in \mathbb{R}^{nxn}$ ortogonales tales que: $A = U\Sigma V^T$ con $\Sigma = diag(\sigma_1, \sigma_2, \dots, \sigma_p) \in \mathbb{R}^{mxn}$, $p = \min\{m,n\}$ y $\sigma_1 \geq \sigma_2 \geq \dots \geq \sigma_p \geq 0$.
Ver 3.3.c.Factorizacion_QR para definición de matriz ortogonal.
Obs: La notación $\sigma_1$ hace referencia al valor singular más grande de A, $\sigma_2$ al segundo valor singular más grande de A y así sucesivamente.
Obs2: La SVD que se definió arriba es nombrada SVD full, hay una forma truncada en la que $U \in \mathbb{R}^{mxk}$, $V \in \mathbb{R}^{nxk}$ y $\Sigma \in \mathbb{R}^{kxk}$.
Existen diferentes propiedades de los valores y vectores singulares, aquí se enlistan algunas:
Si $rank(A) = r$ entonces $r \leq p$ y $\sigma_1 \geq \sigma_2 \geq \dots \geq \sigma_r > \sigma_{r+1} = \sigma_{r+2} = \dots = \sigma_p = 0$.
Si $rank(A) = r$ entonces $A = \displaystyle \sum_{i=0}^r \sigma_i u_i v_i^T$ con $u_i$ $i$-ésima columna de U y $v_i$ $i$-ésima columna de V.
Geométricamente los valores singulares de una matriz $A \in \mathbb{R}^{mxn}$ son las longitudes de los semiejes del hiperelipsoide $E$ definido por $E = \{Ax : ||x|| \leq 1, \text{ con } ||\cdot || \text{ norma Euclidiana}\}$ y los vectores $u_i$ son direcciones de estos semiejes; los vectores $vi$'s tienen norma igual a $1$ por lo que se encuentran en una circunferencia de radio igual a $1$ y como $Av_i = \sigma u_i$ entonces $A$ mapea los vectores $v_i$'s a los semiejes $u_i$'s respectivamente:
Entre las aplicaciones de la SVD se encuentran:
La idea del método de Jacobi one sided consiste en multiplicar a la matriz $A \in \mathbb{R}^{m \times n}$ por la derecha de forma repetida por matrices ortogonales de nombre rotaciones Givens hasta que se converja a $U \Sigma$.
Si $u, v \in \mathbb{R}^2-\{0\}$ con $\ell = ||u||_2 = ||v||_2$ y se desea rotar al vector $u$ en sentido contrario a las manecillas del reloj por un ángulo $\theta$ para llevarlo a la dirección de $v$:
A partir de las relaciones anteriores como $cos(\phi)=\frac{u_1}{\ell}, sen(\phi)=\frac{u_2}{\ell}$ se tiene: $v_1 = (cos\theta)u_1-(sen\theta)u_2$, $v_2=(sen\theta)u_1+(cos\theta)u_2$ equivalentemente:
$$\begin{array}{l} \left[\begin{array}{c} v_1\\ v_2 \end{array} \right] = \left[ \begin{array}{cc} cos\theta & -sen\theta\\ sen\theta & cos\theta \end{array} \right] \cdot \left[\begin{array}{c} u_1\\ u_2 \end{array} \right] \end{array} $$Comentarios:
La matriz $R_O$ se llama matriz de rotación o rotaciones Givens, es una matriz ortogonal pues $R_O^TR_O=I_2$.
La multiplicación $v=R_Ou$ es una rotación en sentido contrario a las manecillas del reloj. La multiplicación $u=R_O^Tv$ es una rotación en sentido de las manecillas del reloj y el ángulo asociado es $-\theta$.
Obsérvese que $det(R_O)=1$.
Ejemplo:
Rotar al vector $v=(1,1)^T$ un ángulo de $45^o$ en sentido contrario a las manecillas del reloj.
Solución:
In [1]:
import numpy as np
import math
In [2]:
v=np.array([1,1])
La matriz $R_O$ es: $R_O = \left[ \begin{array}{cc} cos(\frac{\pi}{4}) & -sen(\frac{\pi}{4})\\ sen(\frac{\pi}{4}) & cos(\frac{\pi}{4}) \end{array} \right ] $
In [3]:
theta=math.pi/4
RO=np.array([[math.cos(theta), -math.sin(theta)],
[math.sin(theta), math.cos(theta)]])
In [8]:
RO
Out[8]:
In [4]:
RO@v
Out[4]:
Obsérvese en el ejemplo anterior que se preservó su la longitud de $v$:
In [5]:
np.linalg.norm(v)
Out[5]:
pues:
$$\begin{array}{l} \left[ \begin{array}{cc} \frac{v_1}{\sqrt{v_1^2+v_2^2}} & \frac{v_2}{\sqrt{v_1^2+v_2^2}}\\ -\frac{v_2}{\sqrt{v_1^2+v_2^2}} & \frac{v_1}{\sqrt{v_1^2+v_2^2}} \end{array} \right ] \cdot \left[\begin{array}{c} v_1\\ v_2 \end{array} \right]= \left[ \begin{array}{c} \frac{v_1^2+v_2^2}{\sqrt{v_1^2+v_2^2}}\\ \frac{-v_1v_2+v_1v_2}{\sqrt{v_1^2+v_2^2}} \end{array} \right ] = \left[ \begin{array}{c} \frac{v_1^2+v_2^2}{\sqrt{v_1^2+v_2^2}}\\ 0 \end{array} \right ]= \left[ \begin{array}{c} ||v||_2\\ 0 \end{array} \right ] \end{array} $$Y definiendo $cos(\theta)=\frac{v_1}{\sqrt{v_1^2+v_2^2}}, sen(\theta)=\frac{v_2}{\sqrt{v_1^2+v_2^2}}$ se tiene :
$$ R_O=\left[ \begin{array}{cc} cos\theta & sen\theta\\ -sen\theta & cos\theta \end{array} \right] $$que en el ejemplo anterior como $v=(1,1)^T$ entonces: $cos(\theta)=\frac{1}{\sqrt{2}}, sen(\theta)=\frac{1}{\sqrt{2}}$ por lo que $\theta=\frac{\pi}{4}$ y:
$$ R_O=\left[ \begin{array}{cc} \frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}}\\ -\frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}} \end{array} \right] $$Para hacer cero la entrada $v_1$ de $v$ hay que usar:
$$\begin{array}{l} R_O=\left[ \begin{array}{cc} cos\theta & -sen\theta\\ sen\theta & cos\theta \end{array} \right] =\left[ \begin{array}{cc} \frac{1}{\sqrt{2}} & -\frac{1}{\sqrt{2}}\\ \frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}} \end{array} \right] \end{array} $$El producto de las rotaciones Givens construye a la matriz ortogonal $V \in \mathbb{R}^{n \times n}$:
$$AV \rightarrow W \in \mathbb{R}^{m \times n}$$Comentarios:
con $U_1 \in \mathbb{R}^{m \times r}$ matriz con columnas ortonormales: $U_1^TU_1=I_r$ y $\Sigma = diag(\sigma_1,\dots, \sigma_r)$ matriz diagonal.
Algoritmo one sided Jacobi
Entrada:
matriz $A \in \mathbb{R}^{m \times n}$: matriz a la que se le calculará su SVD.
$TOL$ controla la convergencia del método.
$maxsweeps$ número máximo de sweeps (descritos en los comentarios de abajo).
Salida:
Nota: se utilizará la notación $A^{(k)}=[a_1^{(k)} a_2^{(k)} \cdots a_n^{(k)}]$ con cada $a_i^{(k)}$ como $i$-ésima columna de $A$ y $k$ representa el sweep.
Definir $A^{(0)}=A$, $V^{(0)}=I_n$ (sweep $=0$).
Mientras no se haya llegado al número máximo de sweeps ($k \leq maxsweeps$) o el número de columnas ortogonales sea menor a $\frac{n(n-1)}{2}$:
Para todos los pares con índices $i<j$ generados con alguna metodología (descrita en la sección de abajo) y para k desde 0 hasta convergencia:
Revisar si las columnas $a_i^{(k)}, a_j^{(k)}$ de $A^{(k)}$ son ortogonales (el chequeo se describe en los comentarios). Si son ortogonales se incrementa por uno la variable $num\text{_}columnas\text{_}ortogonales$. Si no son ortogonales:
Calcular $\left[ \begin{array}{cc} a & c\\ c & b \end{array} \right]$ la submatriz $(i,j)$ de $A^{T(0)}A^{(0)}$ donde: $a = ||a_i^{(k)}||_2^2, b=||a_j^{(k)}||_2^2, c=a_i^{T(k)}a_j^{(k)}$.
Calcular la rotación Givens que diagonaliza $\left[ \begin{array}{cc} a & c\\ c & b \end{array} \right]$ definiendo: $\xi = \frac{b-a}{2c}, t=\frac{signo(\xi)}{|\xi| + \sqrt{1+\xi^2}}, cs = \frac{1}{\sqrt{1+t^2}}, sn = cs*t$. Recordar que $signo(a) = \begin{cases} 1 &\text{ si } a \geq 0 ,\\ -1 &\text{ si } a < 0 \end{cases}$.
Actualizar las columnas $i,j$ de $A^{(k)}$. Para $\ell$ de $1$ a $n$:
$temp = A^{(k)}_{\ell i}$
$A_{\ell i}^{(k)} = cs*temp - sn*A_{\ell j}^{(k)}$
$A_{\ell j}^{(k)} = sn*temp + cs*A_{\ell j}^{(k)}$
Actualizar a la matriz $V^{(k)}$. Para $\ell$ de $1$ a $n$:
$temp = V_{\ell i}^{(k)}$
$V_{\ell i}^{(k)} = cs*temp - sn*V_{\ell j}^{(k)}$
$V_{\ell j}^{(k)} = sn*temp + cs*V_{\ell j}^{(k)}$
Incrementar por uno la variable $k$ que cuenta el número de sweeps.
Comentarios:
Las rotaciones se realizan en una secuencia con nombre sweep. Cada sweep consiste de como máximo $\frac{n(n-1)}{2}$ rotaciones (pues depende de cuántas columnas son o no ortogonales) y en cada sweep se ortogonalizan $2$ columnas. El número de sweeps a realizar se controla con la variable $maxsweeps$ y con la variable $num\text{_}columnas\text{_}ortogonales$ que va contando en cada sweep el número de columnas ortogonales.
La convergencia del algoritmo anterior involucra dos aspectos:
El número de columnas ortogonales que tenemos en un sweep: si en un sweep tal número es $\frac{n(n-1)}{2}$ (almacenado en la variable $num\text{_}columnas\text{_}ortogonales$) entonces el algoritmo termina.
¿Cómo reviso si las columnas $i,j$ de $A^{(k)}$ son ortogonales? si se cumple que
con $TOL$ un valor menor o igual a $10^{-8}$ entonces son ortogonales las columnas $a_i^{(k)}, a_j^{(k)}$ de $A^{(k)}$.
El ángulo $\theta$ se elige de acuerdo a $2$ criterios:
1) El ángulo es $0$ si $a_i^{T(k)}a_j^{(k)}=0$ y por lo tanto las columnas $i,j$ son ortogonales y no se hace rotación.
2) $\theta \in (\frac{-\pi}{4}, \frac{pi}{4})$ tal que $a_i^{T(k+1)}a_j^{(k+1)}=0$ y para este caso se calculan $\xi, t, cs,sn$ (las variables $cs, sn$ hacen referencia a $cos(\theta), sen(\theta)$.
Las actualizaciones para $A,V$ en el algoritmo son de la forma: $A^{(k+1)} = A^{(k)}U^{(k)}, V^{(k+1)} = V^{(k)}U^{(k)}$ para $k>0$ donde: $U^{(k)}$ son matrices de rotación del plano $(i,j)$. Esto es, una identidad pero con elementos:
¿Cómo elijo los pares de columnas (i,j) a las que se les aplicará un proceso de ortogonalización?
Hay distintas metodologías para la elección de las columnas $a_i^{(k)}, a_j^{(k)}$ de $A^{(k)}$. Se presenta a continuación una de ellas:
Paso 1: encontrar factores $U, \Sigma, V$ tales que $A=U \Sigma V^T$.
Paso 2: resolver el sistema diagonal $\Sigma d = U^Tb$.
Paso 3: realizar la multiplicación matriz vector para obtener $x$: $x=Vd$.
Referencias: