!pip install --quiet sklearn

import math

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA

np.random.seed(2020)
mpoints=200
X = (np.random.rand(2,2)@np.random.normal(0,1,(2,mpoints))).T

plt.scatter(X[:, 0], X[:, 1])
plt.grid()
plt.axis('equal');

pca = PCA(n_components=2,svd_solver='full')
pca.fit(X)

"""
Explanation: (VALVECSINGALGSVD)=
2.4 Valores, vectores singulares y algoritmos para calcular la SVD
```{admonition} Notas para contenedor de docker:
Comando de docker para ejecución de la nota de forma local:
nota: cambiar &lt;ruta a mi directorio&gt; por la ruta de directorio que se desea mapear a /datos dentro del contenedor de docker y &lt;versión imagen de docker&gt; por la versión más actualizada que se presenta en la documentación.
docker run --rm -v &lt;ruta a mi directorio&gt;:/datos --name jupyterlab_optimizacion -p 8888:8888 -d palmoreck/jupyterlab_optimizacion:&lt;versión imagen de docker&gt;
password para jupyterlab: qwerty
Detener el contenedor de docker:
docker stop jupyterlab_optimizacion
Documentación de la imagen de docker palmoreck/jupyterlab_optimizacion:&lt;versión imagen de docker&gt; en liga.
```

Nota generada a partir de liga.
```{admonition} Al final de esta nota la comunidad lectora:
:class: tip


Aprenderá algunas definiciones y resultados de los valores y vectores singulares.


Comprenderá la relación entre la descomposición de valores singulares de una matriz y las componentes principales desde una perspectiva algebraica. Así también tendrá definiciones utilizadas en Estadística para este último tema.


Aprenderá la relación entre eigenvalores y las matrices definidas y semidefinidas positivas.


Se proporcionará una lista de algoritmos para calcular la descomposición en valores singulares.


```
En esta nota asumimos que $A \in \mathbb{R}^{m \times n}$.
Valor singular
```{admonition} Definición
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.
```
```{admonition} Observación
:class: tip
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 no negativos 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
```{admonition} Definición
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$.
Por ejemplo para un caso $m < n$:
$$
\left [
\begin{array}{ccc}
a_{11} & a_{12} & a_{13}\
a_{21} & a_{22} & a_{23}
\end{array}
\right ]
=
\left [
\begin{array}{cc}
u_{11} & u_{12}\
u_{21} & u_{22}
\end{array}
\right ]
\left [
\begin{array}{ccc}
\sigma_1 & 0 & 0\
0 & \sigma_2 & 0
\end{array}
\right ]
\left [
\begin{array}{ccc}
v_{11} & v_{21} & v_{31}\
v_{12} & v_{22} & v_{32}\
v_{13} & v_{23} & v_{33}
\end{array}
\right ]
$$
y para un caso $m > n$:
$$
\left [
\begin{array}{cc}
a_{11} & a_{12}\
a_{21} & a_{22}\
a_{31} & a_{32}
\end{array}
\right ]
=
\left [
\begin{array}{ccc}
u_{11} & u_{12} & u_{13}\
u_{21} & u_{22} & u_{23}\
u_{31} & u_{32} & u_{33}
\end{array}
\right ]
\left [
\begin{array}{cc}
\sigma_1 & 0\
0 & \sigma_2\
0 & 0
\end{array}
\right ]
\left [
\begin{array}{cc}
v_{11} & v_{21}\
v_{12} & v_{22}\
\end{array}
\right ]
$$
```{admonition} Definición
Las columnas de $U$ nombramos vectores singulares izquierdos de $A$ y las columnas de $V$ nombramos vectores singulares derechos de $A$ en $A = U \Sigma V^T$.
```
```{admonition} Comentarios


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.


Para cualquier $A \in \mathbb{R}^{m \times n}$ se tiene $A = U \Sigma V^T$. Si $b = Ax = (U \Sigma V^T) x$ entonces:


$$\tilde{b} = U^Tb = U^T  (Ax) = U^T (U \Sigma V^T) x = \Sigma V^Tx = \Sigma \tilde{x}.$$
Lo anterior indica que el producto matricial $Ax$ para cualquier matriz $A$ es equivalente a multiplicar una matriz diagonal por un vector denotado como $\tilde{x}$ que contiene los coeficientes de la combinación lineal de las columnas de $V$ para el vector $x$ . El resultado de tal multiplicación es un vector denotado como $\tilde{b}$ que contiene los coeficientes de la combinación lineal de las columnas de $U$ para el vector $b$. En resúmen, la multiplicación $Ax$ es equivalente a la multiplicación por una matriz diagonal $\Sigma \tilde{x}$ salvo dos cambios de bases, ver Change of basis, la base de los vectores singulares derechos (columnas de $V$) y la base de los vectores singulares izquierdos (columnas de $U$).

La SVD que se definió arriba es nombrada SVD full, hay otras formas como la truncada en la que $U \in \mathbb{R}^{m \times k}$, $V \in \mathbb{R}^{nxk}$ y $\Sigma \in \mathbb{R}^{k \times k}$ con $k \leq r$, compacta donde $k=r$ y la thin en la que $k=p$:

<img src="https://dl.dropboxusercontent.com/s/8dq0jiw5em93k1j/svd_thin.png?dl=0" heigth="700" width="700">
donde: $r = rank(A)$.
Ver reduced SVDs.
```
Algunas propiedades
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||_2 \leq 1}$ 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:


<img src="https://dl.dropboxusercontent.com/s/1yqoe4qibyyej53/svd_2.jpg?dl=0" heigth="700" width="700">

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

<img src="https://dl.dropboxusercontent.com/s/uo9s9f0nqi43s6d/svd_four_spaces_of_matrix.png?dl=0" heigth="600" width="600">

Si $t < r$ y $r=rank(A)$ entonces $A_t =  \displaystyle \sum_{i=0}^t \sigma_i u_i v_i^T$ (SVD truncada) 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 la norma matricial Euclidiana y de Frobenius, en el caso de Frobenius es la única matriz que cumple lo anterior.

Algunas aplicaciones
Algunas de 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.

Ejemplo: componentes principales
Considérese $X \in \mathbb{R}^{m \times n}$ cuyas columnas tienen una observación de un vector aleatorio (tenemos $n$ vectores aleatorios de mediciones) y $X = U \Sigma V^T$ la descomposición en valores singulares de $X$, entonces los vectores singulares derechos $v_i$ (columnas de la matriz $V$) son nombrados ejes o direcciones principales de $X$ y el vector $z_1 = X v_1 = \sigma_1 u_1$ con $u_1$ vector singular izquierdo (primera columna de la matriz $U$) tiene varianza muestral:
$$\text{var}(z_1) = \text{var}(X v_1)= \text{var}(\sigma_1 u_1) = \sigma_1^2 \text{var}(u_1) = \sigma_1^2 \left [ \frac{1}{m} \displaystyle \sum_{i=1}^m (u_1(i) - \bar{u}_1)^2 \right ]$$
donde: $u_1(i)$ es la $i$-ésima componente de $u_1$ y $\sigma_1$ es el máximo valor singular de $X$ también denotado como $\sigma_{\text{max}}$.
```{admonition} Comentarios
Si la media de cada columna de $X$ es cero, $X$ se nombra centrada, entonces:

$z_1$ tiene la máxima varianza muestral entre todas las combinaciones lineales de las columnas de $X$ pues:

$$\text{var}(z_1) = \frac{\sigma_1^2}{m} \displaystyle \sum_{i=1}^m u_1(i)^2 = \frac{\sigma_1^2}{m} ||u_1||_2^2 = \frac{\sigma_1^2}{m}.$$

$z_1$ es la primera componente principal y el vector $u_1 = \frac{1}{\sigma_1}z_1 = \frac{1}{\sigma_1}Xv_1$ se le nombra primera componente principal normalizada. El vector $v_1$ es la primera dirección principal de $X$ o también nombrada loading. 

```
Algunas definiciones utilizadas en Estadística se encuentran las siguientes.
```{admonition} Definición
La matriz de correlaciones entre cada componente principal normalizada $u$'s y cada columna de $X$ es:
$$ C = \left (\frac{\sigma_1}{\sqrt{m}} v_1 \quad \frac{\sigma_2}{\sqrt{m}} v_2 \cdots \quad \frac{\sigma_n}{\sqrt{m}} v_n \right) \in \mathbb{R}^{n \times n}$$
pues si $x_1$ es la primer columna de $X$ entonces:
$$
\begin{eqnarray}
\text{cov}(x_1,u_1) = \text{cov} \left ( \displaystyle \sum_{k=1}^n \sigma_k v_k(1) u_k, u_1 \right ) &=& \displaystyle \sum_{k=1}^n \text{cov} ( \sigma_k v_k(1) u_k, u_1 ) \nonumber \
&=& \displaystyle \sum_{k=1}^n \sigma_k v_k(1) \text{cov} (u_k, u_1) \nonumber \
&=& \sigma_1 v_1(1) \text{var}(u_1) \nonumber \
&=& \frac{\sigma_1 v_1(1)}{m} \sum_{i=1}^m u_1(i)^2 \nonumber \
&=& \frac{\sigma_1 v_1(1)}{m} \nonumber
\end{eqnarray}
$$
Y como $\text{cor}(x_1,u_1) = \frac{\text{cov}(x_1,u_1)}{\sqrt{\text{var}(x_1)} \sqrt{\text{var}(u_1)}}$ se tiene:
$$\text{cor}(x_1,u_1) = \frac{\frac{\sigma_1 v_1(1)}{m}}{1 \cdot \sqrt{\frac{1}{m}}}  = \frac{\sigma_1 v_1(1)}{\sqrt{m}} $$
```
```{admonition} Definición
El cociente de varianza explicada para cada componente es el número:
$$\frac{\sigma_i^2}{\displaystyle \sum_{i=1}^p \sigma_i^2}$$
con $p = \min(m,n)$.
```
```{admonition} Observación
:class: tip
la matriz $\frac{1}{m}X^TX$ es la matriz de varianzas y covarianzas muestral la cual siempre es una matriz simétrica positiva semidefinida (aún si la $X$ no es centrada). Si $X$ además de ser centrada cumple que la varianza de cada una de sus columnas es $1$, $X$ se nombra estandarizada. La matriz $\frac{1}{m}X^TX$ en este caso es la matriz de correlaciones muestral.
```
(MATRIZSDPOSITIVA)=
```{admonition} Definición
Con los eigenvalores podemos caracterizar a las matrices definidas y semidefinidas positivas: $A$ es semidefinida positiva si y sólo si los eigenvalores de $T=\frac{A+A^T}{2}$ son no negativos. Es definida positiva si y sólo si los eigenvalores de $T$ son positivos. 
Los conjuntos de matrices que se utilizan para definir a matrices simétricas semidefinidas positivas y simétricas definidas positivas son $\mathbb{S}{+}^n$ y $\mathbb{S}{++}^n$ respectivamente ($\mathbb{S}$ es el conjunto de matrices simétricas).
```
End of explanation
"""


print(pca.components_)

"""
Explanation: Direcciones o ejes principales:
End of explanation
"""


print(pca.explained_variance_ratio_)

"""
Explanation: Varianza explicada por cada componente (primera posición para la primera componente):
End of explanation
"""


pca.singular_values_

"""
Explanation: Valores singulares de la matriz $X$:
End of explanation
"""


z = pca.transform(X)

z.shape

"""
Explanation: Componentes principales con el método transform de scikit-learn.PCA:
End of explanation
"""


z1 = z[:,0]

z1[0:10]

"""
Explanation: Primera componente:
End of explanation
"""


z2 = z[:,1]

z2[0:10]

"""
Explanation: Segunda componente:
End of explanation
"""


X_centered = X - X.mean(axis=0)

"""
Explanation: Usemos el método de SVD del paquete numpy para revisar lo que nos devuelve y observar que son iguales sus resultados:
Primero centramos a la $X$:
End of explanation
"""


u,s,vt = np.linalg.svd(X_centered)

"""
Explanation: No olvidemos que el método de numpy nos devuelve $V^T$ y no $V$:
End of explanation
"""


print(np.isclose(pca.singular_values_,s))

"""
Explanation: Los valores singulares están dados por la diagonal de $\Sigma$:
End of explanation
"""


print(np.isclose(np.abs(pca.components_[0,:]),np.abs(vt.T[:,0])))

print(np.isclose(np.abs(pca.components_[1,:]), np.abs(vt.T[:,1])))

"""
Explanation: Las direcciones principales están dadas por las columnas de $V$ (salvo signos positivos o negativos):
End of explanation
"""


z_manual = X_centered@(vt.T)

z1_manual = z_manual[:,0]

print(np.isclose(np.abs(z1), np.abs(z1_manual)))

z2_manual = z_manual[:,1]

print(np.isclose(np.abs(z2), np.abs(z2_manual)))

"""
Explanation: Las componentes principales están dadas por la multiplicación matricial $XV$ (salvo signos positivos o negativos):
End of explanation
"""


z1_manual_2 = s[0]*u[:,0]

print(np.isclose(np.abs(z1), np.abs(z1_manual_2)))

z2_manual_2 = s[1]*u[:,1]

print(np.isclose(np.abs(z2), np.abs(z2_manual_2)))

"""
Explanation: También podemos hacer la multiplicación $\sigma u$ (salvo signos positivos o negativos):
End of explanation
"""


print(np.isclose(pca.explained_variance_ratio_[0], s[0]**2/np.sum(s**2)))

print(np.isclose(pca.explained_variance_ratio_[1], s[1]**2/np.sum(s**2)))

"""
Explanation: La varianza explicada está dada por los valores singulares al cuadrado divididos por la suma de éstos al cuadrado:
End of explanation
"""


