<jupyter_start><jupyter_text>(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 `` por la ruta de directorio que se desea mapear a `/datos` dentro del contenedor de docker y `` por la versión más actualizada que se presenta en la documentación.`docker run --rm -v :/datos --name jupyterlab_optimizacion -p 8888:8888 -d palmoreck/jupyterlab_optimizacion:`password para jupyterlab: `<PASSWORD>`Detener el contenedor de docker:`docker stop jupyterlab_optimizacion`Documentación de la imagen de docker `palmoreck/jupyterlab_optimizacion:` en [liga](https://github.com/palmoreck/dockerfiles/tree/master/jupyterlab/optimizacion).``` --- Nota generada a partir de [liga](https://www.dropbox.com/s/s4ch0ww1687pl76/3.2.2.Factorizaciones_matriciales_SVD_Cholesky_QR.pdf?dl=0). ```{admonition} Al final de esta nota el y la 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 definición de matriz definida y semidefinida positiva.* 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ónEl 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: tipLa 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ónAsociado 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ónLas 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](https://en.wikipedia.org/wiki/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$:donde: $r = rank(A)$.Ver [reduced SVDs](https://en.wikipedia.org/wiki/Singular_value_decompositionReduced_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:  * 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$ (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} ComentariosSi 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ónLa **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ónEl 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: tipla 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ónUna matriz $A$ es semidefinida positiva si $x^TAx \geq 0$ $\forall x \in \mathbb{R}^n - \{0\}$. Si se cumple de forma estricta la desigualdad entonces $A$ es **definida positiva**. 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).```<jupyter_code>!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)<jupyter_output><empty_output><jupyter_text>**Direcciones o ejes principales**:<jupyter_code>print(pca.components_)<jupyter_output>[[ 0.91945871  0.39318658]
 [ 0.39318658 -0.91945871]]<jupyter_text>**Varianza explicada por cada componente (primera posición para la primera componente):**<jupyter_code>print(pca.explained_variance_ratio_)<jupyter_output>[0.99029474 0.00970526]<jupyter_text>**Valores singulares de la matriz $X$:**<jupyter_code>pca.singular_values_<jupyter_output><empty_output><jupyter_text>**Componentes principales con el método [transform](https://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.htmlsklearn.decomposition.PCA.transform) de `scikit-learn.PCA`:**<jupyter_code>z = pca.transform(X)
z.shape<jupyter_output><empty_output><jupyter_text>**Primera componente:**<jupyter_code>z1 = z[:,0]
z1[0:10]<jupyter_output><empty_output><jupyter_text>**Segunda componente:**<jupyter_code>z2 = z[:,1]
z2[0:10]<jupyter_output><empty_output><jupyter_text>Usemos el método de [SVD](https://numpy.org/doc/stable/reference/generated/numpy.linalg.svd.htmlnumpy.linalg.svd) del paquete `numpy` para revisar lo que nos devuelve y observar que son iguales sus resultados: Primero centramos a la $X$:<jupyter_code>X_centered = X - X.mean(axis=0)<jupyter_output><empty_output><jupyter_text>No olvidemos que el método de `numpy` nos devuelve $V^T$ y no $V$:<jupyter_code>u,s,vt = np.linalg.svd(X_centered)<jupyter_output><empty_output><jupyter_text>**Los valores singulares están dados por la diagonal de $\Sigma$:**<jupyter_code>print(np.isclose(pca.singular_values_,s))<jupyter_output>[ True  True]<jupyter_text>**Las direcciones principales están dadas por las columnas de $V$ (salvo signos positivos o negativos):**<jupyter_code>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])))<jupyter_output>[ True  True]<jupyter_text>**Las componentes principales están dadas por la multiplicación matricial $XV$ (salvo signos positivos o negativos):**<jupyter_code>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)))<jupyter_output>[ True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  T[...]<jupyter_text>**También podemos hacer la multiplicación $\sigma u$ (salvo signos positivos o negativos):**<jupyter_code>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)))<jupyter_output>[ True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  T[...]<jupyter_text>**La varianza explicada está dada por los valores singulares al cuadrado divididos por la suma de éstos al cuadrado:**<jupyter_code>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)))<jupyter_output>True
{'fraction_non_alphanumeric': 0.13784959602237415, 'fraction_numerical': 0.02144188937228092, 'mean_word_length': 4.65787623066104, 'pattern_counts': {'":': 0, '<': 53, '<?xml version=': 0, '>': 53, 'https://': 6, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '215739', 'license': 'Apache-2.0', 'n_tokens_mistral': 6182, 'n_tokens_neox': 6031, 'n_words': 1879, 'path': 'libro_optimizacion/temas/2.computo_matricial/2.4/Valores_vectores_singulares_y_algoritmos_para_calcular_la_SVD.ipynb', 'repo_name': 'ITAM-DS/analisis-numerico-computo-cientifico', 'size': '16086'}