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

password para jupyterlab: qwerty

Detener el contenedor de docker:

docker stop jupyterlab_numerical

Documentación de la imagen de docker palmoreck/jupyterlab_numerical:1.1.0 en liga.


Nota basada en liga

Definiciones generales

En lo que sigue se supone una matriz cuadrada $A \in \mathbb{R}^{nxn}$.

Eigenvalor (valor propio o característico)

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:

  • Una matriz con componentes reales puede tener eigenvalores y eigenvectores con valores en $\mathbb{C}$ o $\mathbb{C}^n$ respectivamente.
  • El conjunto de eigenvalores se le llama espectro de una matriz.
  • $A$ siempre tiene al menos un eigenvalor con eigenvector asociado.

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

Valores y vectores singulares de una matriz

En lo que sigue se supone $A \in \mathbb{R}^{mxn}$.

Valor singular

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

Vector singular izquierdo, vector singular derecho

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.

Descomposición en valores singulares (SVD)

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:

  • La SVD da bases ortogonales para los $4$ espacios fundamentales de una matriz: espacio columna, espacio nulo izquierdo, espacio nulo y espacio renglón:

  • Si $t < r$ y $r=rank(A)$ entonces $A_t = \displaystyle \sum_{i=0}^t \sigma_i u_i v_i^T$ es una matriz de entre todas las matrices con $rank$ igual a t, que es más cercana a A (la cercanía se mide con una norma matricial).

Entre las aplicaciones de la SVD se encuentran:

  • Procesamiento de imágenes y señales.
  • Sistemas de recomendación (Netflix).
  • Mínimos cuadrados.
  • Componentes principales.
  • Reconstrucción de imágenes.

Método de Jacobi para calcular la SVD

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

Rotaciones

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]:
array([[ 0.70710678, -0.70710678],
       [ 0.70710678,  0.70710678]])

In [4]:
RO@v


Out[4]:
array([1.11022302e-16, 1.41421356e+00])

Obsérvese en el ejemplo anterior que se preservó su la longitud de $v$:


In [5]:
np.linalg.norm(v)


Out[5]:
1.4142135623730951
  • Esto se cumple para todas las matrices ortogonales. Son isometrías bajo la norma $2$ o Euclidiana: $||R_0v||_2=||v||_2$
  • Obsérvese en el ejemplo anterior que se hizo cero la entrada $v_1$ de $v$. Las matrices de rotación se utilizan para hacer ceros en entradas de un vector. Por ejemplo si $v=(v_1,v_2)^T$ y se desea hacer cero la entrada $v_2$ se puede utilizar la matriz de rotación:
$$R_O = \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 ] $$

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} $$

Algoritmo de Jacobi one sided

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:

  • Las normas Euclidianas de las columnas de $W$ construyen a los valores singulares $\sigma_i \forall i=1,\dots,r$:
$$W = [U_1 \quad 0]\left[ \begin{array}{cc} \Sigma & 0\\ 0 & 0 \end{array} \right]$$

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:

  • matrices $V \in \mathbb{R}^{n \times n}$ ortogonal, $W \in \mathbb{R}^{m \times n}$ representada en el algoritmo por $A^{(k)}$ para un valor de $k$ controlado por la convergencia (descrita en los comentarios de abajo).

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

$$\frac{|a_i^{T (k)}a_j^{(k)}|}{||a_i^{(k)}||_2||a_j^{(k)}||_2} < TOL$$

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:

$$u_{ii}^{(k)} = cos(\theta), \quad u_{ij}^{(k)} = sen(\theta)$$
$$u_{ji}^{(k)}=-sen(\theta), \quad u_{jj}^{(k)}=cos(\theta).$$
  • La multiplicación $A^{(k)}U^{(k)}$ sólo involucra a $2$ columnas de $A^{(k)}$:
$$(a_i^{(k+1)} a_j^{(k+1)}) = (a_i^{(k)} a_j^{(k)})\left[ \begin{array}{cc} cos(\theta) & sen(\theta)\\ -sen(\theta) & cos(\theta) \end{array} \right]$$
  • Al finalizar el algoritmo los valores singulares calculados son las normas Euclidianas de cada columna de $A^{(k)}$. Las columnas normalizadas de $A^{(k)}$ son las columnas de $U$.

¿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:

  • Ordenamiento cíclico por renglones: un sweep involucra trabajar en los pares de columnas $(1,2), (1,3), \dots, (1,n), (2,3), (2,4), \dots, (n-1,n)$. Este ordenamiento siempre converge si $|\theta| \leq \frac{\pi}{4}$.

¿Dado el sistema $Ax=b$, $A \in \mathbb{R}^{n \times n}$ cómo se resuelve con la factorización $SVD$?

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:

  • Ver 4_SVD_y_reconstruccion_de_imagenes para definiciones de eigenvalores, eigenvectores y la descomposición en valores singulares o DVS (o SVD por sus siglas en inglés)