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 generada a partir de liga
In [1]:
!pip3 install -q --user scikit-learn
In [2]:
import os
In [3]:
cur_directory = os.getcwd()
In [4]:
dir_alg_python = '/algoritmos/Python'
In [5]:
os.chdir(cur_directory + dir_alg_python)
En esta nota realizamos una descripción de las componentes principales desde la perspectiva de la descomposición en valores singulares de una matrix $X \in \mathbb{R}^{m \times n}$. El objetivo de la nota es la descripción de los problemas de optimización numérica a resolver en componentes principales y el uso de métodos de descenso (ver 4.2.Algoritmos_para_UCO) para resolver tales problemas, en específico el método de Newton (ver 4.2.Metodo_de_Newton_Python). Se comparan los resultados de la clase de Python PCA del paquete scikit-learn con los obtenidos en la implementación hecha por el prof en algoritmos/Python, en específico algoritmos/Python/algorithms_for_uco.py para problemas tipo UCO (Unconstrained Convex Optimization).
Una referencia didáctica para aprender sobre componentes principales la encuentran en making-sense-of-principal-component-analysis-eigenvectors-eigenvalues.
Supóngase que en cada columna de $X \in \mathbb{R}^{m \times n}$ se tiene una observación de un vector aleatorio (tenemos $n$ vectores aleatorios de mediciones) y sea $X = U \Sigma V^T$ la descomposición en valores singulares de $X$ (ver 3.3.d.SVD). 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 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 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}}$.
Comentarios:
Si la media de cada columna de $X$ es cero, $X$ se nombra centrada, entonces:
el cual tiene solución cerrada dada por: $\sigma_1^2 = \displaystyle \max_{v \in \mathbb{R}^n - \{0\}} \frac{v^TX^TXv}{v^Tv}$, $v_1 = \text{argmax}_{v \in \mathbb{R}^n - \{0\}} \frac{v^TX^TXv}{v^Tv}$ (primera columna de $V$).
con $v_1$ la primera dirección principal de $X$. La solución al problema anterior está dada por: $\sigma_2^2 = \displaystyle \max_{v \in \mathbb{R}^n - \{0\}} \frac{v^TX^TXv}{v^Tv}$, $v_2 = \text{argmax}_{v \in \mathbb{R}^n - \{0\}} \frac{v^TX^TXv}{v^Tv}$ (segunda columna de $V$).
pues si $x_1$ es la primer columna de $X$ entonces:
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}} $$con $p = \min(m,n)$.
Obs: 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.
Como se describió en la sección anterior, las componentes principales tienen una relación directa con la descomposición en valores singulares de $X$. Así, los métodos numéricos para su cómputo involucran métodos para calcular la SVD como son:
Algoritmo de Jacobi one sided, ver 3.3.d.SVD y ex-modulo-3-comp-matricial-svd-czammar para una implementación.
Método de la potencia en el que se utiliza el cociente de Rayleigh para acelerar convergencia. Tal cociente toma la forma $\frac{y^TAy}{y^Ty}$ para $y \neq 0$.
Algoritmo QR basada en la factorización QR, ver 3.3.c.Factorizacion_QR y ex-modulo-3-comp-matricial-qr-dapivei para una implementación.
Métodos de descenso aplicados a problemas de optimización. Recuérdese que en los problemas de optimización convexa sin restricciones (Unconstrained Convex Optimization, ver 4.2.Algoritmos_para_UCO) se resuelven problemas de minimización, sin embargo con una modificación sencilla podemos resolver problemas de maximización.
En esta sección utilizamos el método de descenso por la dirección de Newton, ver 4.2.Metodo_de_Newton_Python, para el cálculo de las componentes principales de la matriz $X \in \mathbb{R}^{m \times n}$ y se comparan los resultados con la clase PCA del paquete scikit-learn de Python.
Construyamos una matriz $X \in \mathbb{R}^{m \times n}$ ejemplo y caculemos sus primeras dos componentes principales vía scikit-learn:
In [6]:
import math
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
from utils import compute_error
In [7]:
np.random.seed(2020)
mpoints=200
X = (np.random.rand(2,2)@np.random.normal(0,1,(2,mpoints))).T
In [8]:
X.shape
Out[8]:
In [9]:
plt.scatter(X[:, 0], X[:, 1])
plt.grid()
plt.axis('equal');
In [10]:
pca = PCA(n_components=2,svd_solver='full')
pca.fit(X)
Out[10]:
Direcciones o ejes principales:
In [11]:
print(pca.components_)
Varianza explicada por cada componente (primera posición para la primera componente):
In [12]:
print(pca.explained_variance_ratio_)
Valores singulares de la matriz $X$:
In [13]:
pca.singular_values_
Out[13]:
Componentes principales con el método transform de scikit-learn.PCA:
In [14]:
z = pca.transform(X)
In [15]:
z.shape
Out[15]:
Primera componente:
In [16]:
z1 = z[:,0]
In [17]:
z1[0:10]
Out[17]:
Segunda componente:
In [18]:
z2 = z[:,1]
In [19]:
z2[0:10]
Out[19]:
numpy para revisar lo que nos devuelve y observar que son iguales sus resultados:Primero centramos a la $X$:
In [20]:
X_centered = X - X.mean(axis=0)
No olvidemos que el método de numpy nos devuelve $V^T$ y no $V$:
In [21]:
u,s,vt = np.linalg.svd(X_centered)
Los valores singulares están dados por la diagonal de $\Sigma$:
In [22]:
compute_error(pca.singular_values_,s)
Out[22]:
Las direcciones principales están dadas por las columnas de $V$ (salvo signos positivos o negativos):
In [23]:
compute_error(np.abs(pca.components_[0,:]),np.abs(vt.T[:,0]), )
Out[23]:
In [24]:
compute_error(np.abs(pca.components_[1,:]),np.abs(vt.T[:,1]))
Out[24]:
Las componentes principales están dadas por la multiplicación matricial $XV$ (salvo signos positivos o negativos):
In [25]:
z_manual = X_centered@(vt.T)
In [26]:
z1_manual = z_manual[:,0]
In [27]:
compute_error(np.abs(z1), np.abs(z1_manual))
Out[27]:
In [28]:
z2_manual = z_manual[:,1]
In [29]:
compute_error(np.abs(z2), np.abs(z2_manual))
Out[29]:
También podemos hacer la multiplicación $\sigma u$ (salvo signos positivos o negativos):
In [30]:
z1_manual_2 = s[0]*u[:,0]
In [31]:
compute_error(np.abs(z1), np.abs(z1_manual_2))
Out[31]:
In [32]:
z2_manual_2 = s[1]*u[:,1]
In [33]:
compute_error(np.abs(z2), np.abs(z2_manual_2))
Out[33]:
La varianza explicada está dada por los valores singulares al cuadrado divididos por la suma de éstos al cuadrado:
In [34]:
compute_error(pca.explained_variance_ratio_[0], s[0]**2/np.sum(s**2))
Out[34]:
In [35]:
compute_error(pca.explained_variance_ratio_[1], s[1]**2/np.sum(s**2))
Out[35]:
El problema de optimización como se planteó al inicio es de la forma:
$$\displaystyle \max_{v \in \mathbb{R}^n - \{0\}} \frac{v^TX^TXv}{v^Tv}$$el cual es equivalente a:
pues se puede asociar $\frac{Xv}{v}$ como $X\hat{v}$ y definir la restricción igual a $1$ anterior.
Como primer enfoque aproximamos a $v_p$ y $\sigma_p$ con $p=\min(m,n)$.
La función objetivo para el mínimo valor singular elevado al cuadrado, $\sigma_{\text{min}}^2$, es:
$$\displaystyle \min_{v \in \mathbb{R}^n - \{0\}} \frac{v^TX^TXv}{v^Tv}$$que es equivalente a:
$$\displaystyle \min_{v \in \mathbb{R}^n} v^TX^TXv$$cuya solución es $\lambda_\text{min}(X^TX)$ y representa el mínimo eigenvalor de $X^TX$ que es igual a $\sigma_{\text{min}}^2$ (mínimo valor singular de $X$ elevado al cuadrado).
Obs: obsérvese que el problema anterior es un problema con función convexa y con restricciones (ver el apéndice de definiciones para funciones convexas de 4.1.Optimizacion_numerica_y_machine_learning). Sin embargo puede reformularse para preservar la convexidad y no tener restricciones, de esta forma utilizaremos algoritmos para problemas convexos sin restricciones (Unconstrained Convex Optimization) revisados en 4.2.Algoritmos_para_UCO. El problema convexo sin restricciones es:
El problema anterior es equivalente al planteado al inicio pues si $f_o(v, \lambda)$ es la función objetivo, entonces:
$$\nabla f_o(v,\lambda) = \left[ \begin{array}{c} X^TX v - \lambda v \\ -\frac{1}{2}(v^Tv-1) \end{array} \right] $$y si igualamos a cero esta ecuación se tiene:
$$\nabla f_o(v,\lambda) = \left[ \begin{array}{c} X^TX v - \lambda v \\ -\frac{1}{2}(v^Tv-1) \end{array} \right] = 0. $$Obs: obsérvese que la variable de optimización es el vector $(v, \lambda) \in \mathbb{R}^{n+1}$.
La ecuación anterior implica: $X^TXv = \lambda v$ y también $v^Tv = 1$ por lo que $(v, \lambda)$ es una pareja (eigenvector, eigenvalor) de la matriz $X^TX$, de hecho el eigenvector tiene norma $1$. Como queremos minimizar $f_o(v,\lambda)$ entonces el óptimo será $\lambda_\text{min}(X^TX)$ con $v_{\lambda\text{min}}$ eigenvector asociado. Así, aproximaremos a $\sigma_{\text{min}}^2$.
Comentario: La formulación anterior tiene una matriz Hessiana dada por la expresión:
$$\nabla^2 f_o(v,\lambda) = \left [ \begin{array}{cc} X^TX - \lambda I & -v \\ -v^T & 0 \end{array} \right ] $$la cual involucra a la matriz $X^TX$. Por este motivo las aproximaciones numéricas que realicemos con este enfoque serán altamente sensibles al número de condición de $X$ y por tanto de $X^TX$.
In [36]:
from algorithms_for_uco import Newtons_method
Función objetivo (que considera a la matriz $X$ centrada):
In [37]:
def fo(x):
v = x[0:(x.size-1)]
value = x[x.size-1]
matvec = X_centered@v
return 1/2*(matvec.dot(matvec)-value*(v.dot(v)-1))
Obs: También podríamos haber definido a la función objetivo con el siguiente código:
def fo(x):
v = x[0:(x.size-1)]
value = x[x.size-1]
m,n=X_centered.shape
cov = (m-1)*np.cov(X_centered, rowvar=False)
return 1/2*(v.dot(cov@v)-value*(v.dot(v)-1))donde estamos usando a la matriz de varianzas y covarianzas dada por la variable cov.
In [38]:
x_ast=np.concatenate((pca.components_[:,1],np.array([pca.singular_values_[1]**2])))
$x^*$:
In [39]:
x_ast
Out[39]:
Punto inicial $x^{(0)}$:
In [40]:
x_0 = np.array([1,0,5],dtype=float)
In [41]:
x_0
Out[41]:
Valor óptimo $p^*$:
In [42]:
p_ast=fo(x_ast)
In [43]:
p_ast
Out[43]:
Argumentos para el método de Newton:
In [44]:
tol=1e-8
tol_backtracking=1e-14
maxiter=50
In [45]:
[x,total_of_iterations,Err_plot,x_plot]=Newtons_method(fo, x_0,tol,
tol_backtracking, x_ast, p_ast, maxiter)
$x$ aproximada por el método de Newton:
In [46]:
x
Out[46]:
$x^*$:
In [47]:
x_ast
Out[47]:
El valor singular $\sigma_{\text{min}}$ es la raíz cuadrada del eigenvalor mínimo de $X^TX$, $\sqrt{\lambda_{\text{min}}(X^TX)}$, que acabamos de calcular y para este ejemplo está en la entrada $3$ del array $x$:
In [48]:
np.sqrt(x[2])
Out[48]:
Error relativo:
In [49]:
compute_error(np.abs(x_ast[0:(x_ast.size-1)]),np.abs(x[0:(x.size-1)]))
Out[49]:
Tenemos alrededor de $6$ dígitos de precisión a la dirección principal $v_{\text{min}}$.
In [50]:
compute_error(x_ast[x_ast.size-1],x[x.size-1])
Out[50]:
Tenemos alrededor de $5$ dígitos de precisión al valor singular $\sigma_{\text{min}}$.
Para la componente principal $z_{\text{min}}$ calculamos el error con:
In [51]:
z = pca.transform(X)
In [52]:
z_approx = X_centered@(x[0:(x.size-1)])
In [53]:
compute_error(np.abs(z[:,1]), np.abs(z_approx))
Out[53]:
Tenemos alrededor de $6$ dígitos de precisión el cálculo de la primera componente principal $z_{\text{min}}$.
Comentario: el método anterior es altamente dependiente al punto inicial, $x^{(0)}$, compárese con otros puntos iniciales y obsérvese este efecto.
Por ejemplo para el punto inicial utilizamos la siguiente definición basada en el cociente de Rayleigh:
Punto inicial $x^{(0)}$:
In [54]:
vec = np.array([0,1],dtype=float) #vector de norma 1
In [55]:
vec
Out[55]:
In [56]:
mv = X_centered@vec
Cociente de Rayleigh:
In [57]:
l = mv.dot(mv)/vec.dot(vec) #para este caso como el vector es de norma 1
# es equivalente hacer: l = mv.dot(mv)
In [58]:
x_0 = np.concatenate((vec,np.array([l])))
In [59]:
x_0
Out[59]:
Valor óptimo $p^*$:
In [60]:
p_ast=fo(x_ast)
In [61]:
p_ast
Out[61]:
Argumentos para el método de Newton:
In [62]:
tol=1e-8
tol_backtracking=1e-14
maxiter=50
In [63]:
[x,total_of_iterations,Err_plot,x_plot] = Newtons_method(fo, x_0,tol,
tol_backtracking, x_ast, p_ast, maxiter)
$x$ aproximada por el método de Newton:
In [64]:
x
Out[64]:
$x^*$:
In [65]:
x_ast
Out[65]:
obsérvese los signos alternados.
Error relativo:
In [66]:
compute_error(np.abs(x_ast[0:(x_ast.size-1)]),np.abs(x[0:(x.size-1)]))
Out[66]:
Tenemos alrededor de $3$ dígitos de precisión a la dirección principal $v_1$.
In [67]:
compute_error(x_ast[x_ast.size-1],x[x.size-1])
Out[67]:
Tenemos alrededor de $2$ dígitos de precisión al valor singular $\sigma_{\text{min}}$.
Para la componente principal $z_{\text{min}}$ calculamos el error con:
In [68]:
z_approx = X_centered@(x[0:(x.size-1)])
In [69]:
compute_error(np.abs(z[:,1]), np.abs(z_approx))
Out[69]:
Tenemos alrededor de $3$ dígitos de precisión el cálculo de la componente principal $z_{\text{min}}$.
Otro Punto inicial $x^{(0)}$:
In [70]:
vec = np.array([10,10],dtype=float) #vector de norma distinta a 1
In [71]:
vec
Out[71]:
In [72]:
mv = X_centered@vec
Definición por cociente de Rayleigh:
In [73]:
l = mv.dot(mv)/vec.dot(vec)
In [74]:
x_0 = np.concatenate((vec,np.array([l])))
In [75]:
x_0
Out[75]:
Valor óptimo $p^*$:
In [76]:
p_ast=fo(x_ast)
In [77]:
p_ast
Out[77]:
Argumentos para el método de Newton:
In [78]:
tol=1e-8
tol_backtracking=1e-14
maxiter=50
In [79]:
[x,total_of_iterations,Err_plot,x_plot] = Newtons_method(fo, x_0,tol,
tol_backtracking, x_ast, p_ast, maxiter)
$x$ aproximada por el método de Newton:
In [80]:
x
Out[80]:
$x^*$:
In [81]:
x_ast
Out[81]:
Error relativo:
In [82]:
compute_error(np.abs(x_ast[0:(x_ast.size-1)]),np.abs(x[0:(x.size-1)]))
Out[82]:
Tenemos $1083\%$ de error a la dirección principal $v_{\text{min}}$.
In [83]:
compute_error(x_ast[x_ast.size-1],x[x.size-1])
Out[83]:
Tenemos $9056\%$ de error al valor singular $\sigma_{\text{min}}$.
Para la componente principal $z_{\text{min}}$ calculamos el error con:
In [84]:
z_approx = X_centered@(x[0:(x.size-1)])
In [85]:
compute_error(np.abs(z[:,1]), np.abs(z_approx))
Out[85]:
Tenemos $11153\%$ de error en el cálculo de la componente principal $z_{\text{min}}$.
La función objetivo para el máximo valor singular elevado al cuadrado, $\sigma_{\text{max}}^2$ es:
$$\displaystyle \max_{v \in \mathbb{R}^n - \{0\}} \frac{v^TX^TXv}{v^Tv}$$Por el desarrollo anterior y la relación entre problemas de minimización y maximización es equivalente resolver:
$$\displaystyle \min_{v \in \mathbb{R}^n, \lambda \in \mathbb{R}} - \left ( \frac{1}{2}(v^TX^TXv - \lambda(v^Tv-1)) \right )$$Función objetivo:
In [86]:
def fo_max(x):
return -fo(x)
In [87]:
x_ast=np.concatenate((pca.components_[:,0],np.array([pca.singular_values_[0]**2])))
$x^*$:
In [88]:
x_ast
Out[88]:
In [89]:
p_ast=fo_max(x_ast)
Valor óptimo $p^*$:
In [90]:
p_ast
Out[90]:
Punto inicial $x^{(0)}$:
In [91]:
x_0 = np.array([1,1,500], dtype=float)
In [92]:
x_0
Out[92]:
Argumentos para el método de Newton:
In [93]:
tol=1e-8
tol_backtracking=1e-14
maxiter=50
In [94]:
[x,total_of_iterations,Err_plot,x_plot] = Newtons_method(fo_max, x_0,tol,
tol_backtracking, x_ast, p_ast, maxiter)
$x$ aproximada por el método de Newton:
In [95]:
x
Out[95]:
$x^*$:
In [96]:
x_ast
Out[96]:
El valor singular $\sigma_{\text{max}}$ es la raíz cuadrada del eigenvalor máximo de $X^TX$, $\sqrt{\lambda_{\text{max}}(X^TX)}$, que acabamos de calcular y para este ejemplo está en la entrada $3$ del array $x$:
In [97]:
np.sqrt(x[2])
Out[97]:
Error relativo:
In [98]:
compute_error(np.abs(x_ast[0:(x_ast.size-1)]),np.abs(x[0:(x.size-1)]))
Out[98]:
Tenemos alrededor de $3$ dígitos de precisión a la dirección principal $v_1$.
In [99]:
compute_error(x_ast[x_ast.size-1],x[x.size-1])
Out[99]:
Tenemos alrededor de $5$ dígitos de precisión al valor singular $\sigma_1$.
Para la componente principal $z_1$ calculamos el error con:
In [100]:
z_approx = X_centered@(x[0:(x.size-1)])
In [101]:
compute_error(np.abs(z[:,0]), np.abs(z_approx))
Out[101]:
Tenemos alrededor de $3$ dígitos de precisión el cálculo de la primera componente principal $z_1$.
Comentario: al igual que en el caso del mínimo valor singular, el método anterior es altamente dependiente al punto inicial, $x^{(0)}$, compárese con otros puntos iniciales y obsérvese este efecto.
Otro Punto inicial $x^{(0)}$:
In [102]:
vec = np.array([1,1],dtype=float) #vector de norma distinta a 1
In [103]:
vec
Out[103]:
In [104]:
mv = X_centered@vec
Definición por cociente de Rayleigh:
In [105]:
l = mv.dot(mv)/vec.dot(vec)
In [106]:
x_0 = np.concatenate((vec,np.array([l])))
In [107]:
x_0
Out[107]:
Valor óptimo $p^*$:
In [108]:
p_ast=fo(x_ast)
In [109]:
p_ast
Out[109]:
Argumentos para el método de Newton:
In [110]:
tol=1e-8
tol_backtracking=1e-14
maxiter=50
In [111]:
[x,total_of_iterations,Err_plot,x_plot] = Newtons_method(fo, x_0,tol,
tol_backtracking, x_ast, p_ast, maxiter)
$x$ aproximada por el método de Newton:
In [112]:
x
Out[112]:
$x^*$:
In [113]:
x_ast
Out[113]:
Error relativo:
In [114]:
compute_error(np.abs(x_ast[0:(x_ast.size-1)]),np.abs(x[0:(x.size-1)]))
Out[114]:
Tenemos $61\%$ de error a la dirección principal $v_1$.
In [115]:
compute_error(x_ast[x_ast.size-1],x[x.size-1])
Out[115]:
Tenemos $13\%$ de error al valor singular $\sigma_1$.
Para la componente principal $z_1$ calculamos el error con:
In [116]:
z_approx = X_centered@(x[0:(x.size-1)])
In [117]:
compute_error(np.abs(z[:,0]), np.abs(z_approx))
Out[117]:
Tenemos $31\%$ de error para el cálculo de la componente principal $z_1$.
Hasta este punto propusimos una matriz con $2$ columnas lo cual ha servido para ilustrar el problema de optimización a resolver y el desempeño del método de Newton para resolverlo. Se vio que la convergencia es muy dependiente del punto inicial y la solución del sistema en el método de Newton que involucra a la Hessiana depende del número de condición de la misma. En esta sección usamos una matriz de más de $2$ columnas y utilizamos un método de deflation clásico para aproximar la segunda componente principal $z_2$.
El método de deflation que utilizamos para aproximar la segunda componente principal $z_2$ consiste en aproximar $z_1$ y $\sigma_1$ como en la sección anterior y posteriormente resolver el problema:
donde: $X_2 = X_1 - \sigma_1 u_1 v_1^T$ y $X_1 = X$.
La $i+1$-ésima componente principal $z_{i+1}$ consiste en resolver el problema:
$$\displaystyle \max_{v \in \mathbb{R}^n - \{0\}} \frac{v^TX_{i+1}^TX_{i+1}v}{v^Tv}$$donde: $X_{i+1} = X_i - \sigma_i u_i v_i^T$ $\forall i=2,3,\dots, p$ con $p = \min(m,n)$.
Tales problemas se resuelven al reescribirlos como:
$$\displaystyle \min_{v \in \mathbb{R}^n, \lambda \in \mathbb{R}} - \left ( \frac{1}{2}(v^TX_i^TX_iv - \lambda(v^Tv-1)) \right )$$para $i=1,2,\dots,p$.
In [118]:
np.random.seed(2020)
mpoints=200
ncols = 3
X = (np.random.rand(ncols,ncols)@np.random.normal(0,1,(ncols,mpoints))).T
In [119]:
pca = PCA(n_components=2,svd_solver='full')
pca.fit(X)
Out[119]:
In [120]:
X_centered = X - X.mean(axis=0)
In [121]:
x_ast=np.concatenate((pca.components_[0,:],np.array([pca.singular_values_[0]**2])))
$x^*:$
In [122]:
x_ast
Out[122]:
Valor óptimo $p^*:$
In [123]:
p_ast=fo_max(x_ast)
In [124]:
p_ast
Out[124]:
Punto inicial $x^{(0)}$:
In [125]:
x_0 = np.array([1,0,1,400], dtype=float)
In [126]:
tol=1e-8
tol_backtracking=1e-14
maxiter=50
In [127]:
[x,total_of_iterations,Err_plot,x_plot] = Newtons_method(fo_max, x_0,tol,
tol_backtracking, x_ast, p_ast, maxiter)
$x$ aproximada por el método de Newton:
In [128]:
x
Out[128]:
$x^*$:
In [129]:
x_ast
Out[129]:
Error relativo:
In [130]:
compute_error(np.abs(x_ast[0:(x_ast.size-1)]),np.abs(x[0:(x.size-1)]))
Out[130]:
Tenemos alrededor de $2$ dígitos de precisión a la dirección principal $v_1$.
In [131]:
compute_error(x_ast[x_ast.size-1],x[x.size-1])
Out[131]:
Tenemos alrededor de $2$ dígitos de precisión al valor singular $\sigma_{1}$.
Para la componente principal $z_1$ calculamos el error con:
In [132]:
z = pca.transform(X)
In [133]:
z_approx = X_centered@(x[0:(x.size-1)])
In [134]:
compute_error(np.abs(z[:,0]), np.abs(z_approx))
Out[134]:
Tenemos alrededor de $2$ dígitos de precisión para la primera componente principal $z_1$.
In [135]:
X_def = X-np.outer(z_approx,x[0:(x.size-1)])
Centramos:
In [136]:
X_centered = X_def - X_def.mean(axis=0)
In [137]:
x_ast=np.concatenate((pca.components_[1,:],np.array([pca.singular_values_[1]**2])))
$x^*$:
In [138]:
x_ast
Out[138]:
Valor óptimo $p^*$:
In [139]:
p_ast
Out[139]:
Punto inicial $x^{(0)}$:
In [140]:
x_0 = np.array([-1,0,0,60], dtype=float)
In [141]:
x_0
Out[141]:
In [142]:
tol=1e-8
tol_backtracking=1e-14
maxiter=50
p_ast=fo_max(x_ast)
In [143]:
[x,total_of_iterations,Err_plot,x_plot] = Newtons_method(fo_max, x_0,tol,
tol_backtracking, x_ast, p_ast, maxiter)
$x$ aproximada por el método de Newton:
In [144]:
x
Out[144]:
$x^*$:
In [145]:
x_ast
Out[145]:
Error relativo:
In [146]:
compute_error(np.abs(x_ast[0:(x_ast.size-1)]),np.abs(x[0:(x.size-1)]))
Out[146]:
Tenemos alrededor de $3$ dígitos de precisión a la dirección principal $v_2$.
In [147]:
compute_error(x_ast[x_ast.size-1],x[x.size-1])
Out[147]:
Tenemos alrededor de $5$ dígitos de precisión al valor singular $\sigma_{2}$.
Para la componente principal $z_2$ calculamos el error con:
In [148]:
z_approx = X_centered@(x[0:(x.size-1)])
In [149]:
compute_error(np.abs(z[:,1]), np.abs(z_approx))
Out[149]:
Tenemos alrededor de $2$ dígitos de precisión para la segunda componente principal $z_2$.
Para este caso utilizamos el gradiente y la Hessiana de forma simbólica en lugar de las aproximaciones por diferencias finitas a la primera y segunda derivada.
Recordamos que las expresiones para el gradiente y la Hessiana para un problema de minimización están dadas por:
Comentario: La formulación anterior tiene una matriz Hessiana dada por la expresión:
$$\nabla^2 f_o(v,\lambda) = \left [ \begin{array}{cc} X^TX - \lambda I & -v \\ -v^T & 0 \end{array} \right ] $$la cual involucra a la matriz $X^TX$. Por este motivo las aproximaciones numéricas que realicemos con este enfoque serán altamente sensibles al número de condición de $X$ y por tanto de $X^TX$.
Gradiente:
In [150]:
def gfo(x):
v = x[0:(x.size-1)]
value = x[x.size-1]
m,n=X_centered.shape
first_block = cov@v-value*v
second_block = -1/2*(v.dot(v)-1)
return np.concatenate((first_block,np.array([second_block])))
In [151]:
def gfo_max(x):
return -gfo(x)
Hessiana:
In [152]:
def Hfo(x):
v = x[0:(x.size-1)]
value = x[x.size-1]
m,n=X_centered.shape
first_block = cov - value*np.eye(v.size)
second_block = -v
fs_block = np.column_stack((first_block, second_block))
third_block = -v.T
fourth_block = np.zeros(1)
tf_block = np.row_stack((third_block.reshape(1,v.size).T, fourth_block)).reshape(1,v.size+1)[0]
return np.row_stack((fs_block, tf_block))
In [153]:
def Hfo_max(x):
return -Hfo(x)
In [154]:
np.random.seed(1)
mpoints=100
ncols=20
X = (np.random.rand(ncols,ncols)@np.random.normal(0,1,(ncols,mpoints))).T
In [155]:
pca = PCA(n_components=2,svd_solver='full')
pca.fit(X)
Out[155]:
In [156]:
X_centered = X - X.mean(axis=0)
In [157]:
m,n = X.shape
Matriz de covarianzas:
In [158]:
cov = (m-1)*np.cov(X_centered, rowvar=False)
In [159]:
x_ast=np.concatenate((pca.components_[0,:],np.array([pca.singular_values_[0]**2])))
$x^*:$
In [160]:
x_ast[0:10]
Out[160]:
In [161]:
x_ast.shape
Out[161]:
In [162]:
x_ast[x_ast.size-1]
Out[162]:
Valor óptimo $p^*:$
In [163]:
p_ast=fo_max(x_ast)
In [164]:
p_ast
Out[164]:
Punto inicial $x^{(0)}$:
In [165]:
x_0 = np.zeros(ncols+1)
In [166]:
x_0[0]=1
In [167]:
x_0[x_0.size-1] = 11000
In [168]:
tol=1e-8
tol_backtracking=1e-14
maxiter=50
In [169]:
[x,total_of_iterations,Err_plot,x_plot] = Newtons_method(fo_max, x_0,tol,
tol_backtracking, x_ast, p_ast, maxiter,
gfo_max,
Hfo_max)
$x$ aproximada por el método de Newton:
In [170]:
x[0:10]
Out[170]:
$x^*$:
In [171]:
x_ast[0:10]
Out[171]:
Error relativo:
In [172]:
compute_error(np.abs(x_ast[0:(x_ast.size-1)]),np.abs(x[0:(x.size-1)]))
Out[172]:
Tenemos alrededor de $6$ dígitos de precisión a la dirección principal $v_1$.
In [173]:
compute_error(x_ast[x_ast.size-1],x[x.size-1])
Out[173]:
Tenemos alrededor de $7$ dígitos de precisión al valor singular $\sigma_1$.
Para la componente principal $z_1$ calculamos el error con:
In [174]:
z = pca.transform(X)
In [175]:
z_approx = X_centered@(x[0:(x.size-1)])
In [176]:
compute_error(np.abs(z[:,0]), np.abs(z_approx))
Out[176]:
Tenemos alrededor de $6$ dígitos de precisión para la primera componente principal $z_1$.
In [177]:
X_def = X-np.outer(z_approx,x[0:(x.size-1)])
Centramos:
In [178]:
X_centered = X_def - X_def.mean(axis=0)
Matriz de covarianzas:
In [179]:
cov = (m-1)*np.cov(X_centered, rowvar=False)
In [180]:
x_ast=np.concatenate((pca.components_[1,:],np.array([pca.singular_values_[1]**2])))
$x^*$:
In [181]:
x_ast[0:10]
Out[181]:
Valor óptimo $p^*$:
In [182]:
p_ast
Out[182]:
In [183]:
x_ast[x_ast.size-1]
Out[183]:
Punto inicial $x^{(0)}$:
In [184]:
x_0 = np.zeros(ncols+1)
In [185]:
x_0[vec.size-1]=1
In [186]:
x_0[x_0.size-1] = 767
In [187]:
tol=1e-8
tol_backtracking=1e-14
maxiter=50
p_ast=fo_max(x_ast)
In [188]:
[x,total_of_iterations,Err_plot,x_plot] = Newtons_method(fo_max, x_0,tol,
tol_backtracking, x_ast, p_ast, maxiter,
gfo_max,
Hfo_max)
$x$ aproximada por el método de Newton:
In [189]:
x[0:10]
Out[189]:
$x^*$:
In [190]:
x_ast[0:10]
Out[190]:
Error relativo:
In [191]:
compute_error(np.abs(x_ast[0:(x_ast.size-1)]),np.abs(x[0:(x.size-1)]))
Out[191]:
Tenemos alrededor de $7$ dígitos de precisión a la dirección principal $v_2$.
In [192]:
compute_error(x_ast[x_ast.size-1],x[x.size-1])
Out[192]:
Tenemos alrededor de $8$ dígitos de precisión al valor singular $\sigma_{2}$.
Para la componente principal $z_2$ calculamos el error con:
In [193]:
z_approx = X_centered@(x[0:(x.size-1)])
In [194]:
compute_error(np.abs(z[:,1]), np.abs(z_approx))
Out[194]:
Tenemos alrededor de $7$ dígitos de precisión para la segunda componente principal $z_2$.
En esta sección calculamos las componentes principales con el método de Newton basándonos en el trabajo realizado por Irene Ramos sobre manejo agrícola con datos de INEGI para municipios de México (ver la liga de su github para el notebook con su trabajo)
In [195]:
# leer los datos
datos = pd.read_csv(cur_directory+ '/datos/censo.csv', encoding ='ISO-8859-1')
datos = datos.set_index('Clave')
In [196]:
datos[['Entidad y municipio', 'Entidad', '% Temporal', '% Mecánica',
'% labor agricola', '% bosque', '% Herbicidas químicos',
'% Insecticidas químicos']].head()
Out[196]:
In [197]:
pca = PCA(n_components=2,svd_solver='full')
In [198]:
datos_subset = datos[['% Temporal', '% Mecánica',
'% labor agricola', '% bosque',
'% Herbicidas químicos',
'% Insecticidas químicos']]
In [199]:
datos_subset.head()
Out[199]:
In [200]:
pca.fit(datos_subset)
Out[200]:
In [201]:
pca.components_
Out[201]:
In [202]:
X = datos_subset.to_numpy()
Centramos:
In [203]:
X_centered = (X - X.mean(axis=0))
In [204]:
m,n = X.shape
Matriz de covarianzas:
In [205]:
cov = (m-1)*np.cov(X_centered, rowvar=False)
$x^*$:
In [206]:
x_ast=np.concatenate((pca.components_[0,:],np.array([pca.singular_values_[0]**2])))
In [207]:
x_ast
Out[207]:
$p^*$:
In [208]:
p_ast=fo_max(x_ast)
In [209]:
p_ast
Out[209]:
Punto inicial $x^{(0)}$:
In [210]:
x_0 = np.zeros(n+1)
In [211]:
x_0[0]=1
In [212]:
x_0[x_0.size-1]=3e6
In [213]:
tol=1e-8
tol_backtracking=1e-14
maxiter=50
In [214]:
[x,total_of_iterations,Err_plot,x_plot] = Newtons_method(fo_max, x_0,tol,
tol_backtracking, x_ast, p_ast, maxiter,
gfo_max,
Hfo_max)
$x$ aproximada por el método de Newton:
In [215]:
x[0:10]
Out[215]:
$x^*$:
In [216]:
x_ast[0:10]
Out[216]:
Error relativo:
In [217]:
compute_error(np.abs(x_ast[0:(x_ast.size-1)]),np.abs(x[0:(x.size-1)]))
Out[217]:
Tenemos alrededor de $2$ dígitos de precisión a la dirección principal $v_1$.
In [218]:
compute_error(x_ast[x_ast.size-1],x[x.size-1])
Out[218]:
Tenemos alrededor de $3$ dígitos de precisión al valor singular $\sigma_{1}$.
Se observa que tienen signos intercambiados $x^*$ y $x$. Se multiplica por $-1$ a $x$:
In [219]:
x = -x
Para la componente principal $z_1$ calculamos el error con:
In [220]:
z = pca.transform(X)
In [221]:
z_approx_1 = X_centered@(x[0:(x.size-1)])
In [222]:
compute_error(np.abs(z[:,0]), np.abs(z_approx_1))
Out[222]:
Tenemos alrededor de $2$ dígitos de precisión para la primera componente principal $z_1$.
In [223]:
X_def = X-np.outer(z_approx_1,x[0:(x.size-1)])
Centramos:
In [224]:
X_centered = X_def - X_def.mean(axis=0)
Matriz de covarianzas:
In [225]:
cov = (m-1)*np.cov(X_centered, rowvar=False)
In [226]:
x_ast=np.concatenate((pca.components_[1,:],np.array([pca.singular_values_[1]**2])))
$x^*$:
In [227]:
x_ast[0:10]
Out[227]:
Valor óptimo $p^*$:
In [228]:
p_ast
Out[228]:
In [229]:
x_ast[x_ast.size-1]
Out[229]:
Punto inicial $x^{(0)}$:
In [230]:
x_0 = np.zeros(n+1)
In [231]:
x_0[x_0.size-2]=1
In [232]:
x_0[x_0.size-1] = 1600000
In [233]:
tol=1e-8
tol_backtracking=1e-14
maxiter=50
p_ast=fo_max(x_ast)
In [234]:
[x,total_of_iterations,Err_plot,x_plot] = Newtons_method(fo_max, x_0,tol,
tol_backtracking, x_ast, p_ast, maxiter,
gfo_max,
Hfo_max)
$x$ aproximada por el método de Newton:
In [235]:
x[0:10]
Out[235]:
$x^*$:
In [236]:
x_ast[0:10]
Out[236]:
Error relativo:
In [237]:
compute_error(np.abs(x_ast[0:(x_ast.size-1)]),np.abs(x[0:(x.size-1)]))
Out[237]:
Tenemos alrededor de $6$ dígitos de precisión a la dirección principal $v_2$.
In [238]:
compute_error(x_ast[x_ast.size-1],x[x.size-1])
Out[238]:
Tenemos alrededor de $8$ dígitos de precisión al valor singular $\sigma_{2}$.
Para la componente principal $z_2$ calculamos el error con:
In [239]:
z_approx_2 = X_centered@(x[0:(x.size-1)])
In [240]:
compute_error(np.abs(z[:,1]), np.abs(z_approx_2))
Out[240]:
Tenemos alrededor de $6$ dígitos de precisión a la componente principal $z_2$.
In [241]:
from sklearn.cluster import KMeans
import matplotlib.colors
In [242]:
z_approx = np.column_stack((z_approx_1, z_approx_2))
In [243]:
z_approx #estas son las dos componentes principales
Out[243]:
Método de K-means
In [244]:
kmeans = KMeans(n_clusters = 4)
kmeans.fit(z_approx)
centers = kmeans.cluster_centers_
labels = kmeans.labels_
In [245]:
gold = matplotlib.colors.to_rgba('Gold')
green = matplotlib.colors.to_rgba('SeaGreen')
purple = matplotlib.colors.to_rgba('Indigo')
blue = matplotlib.colors.to_rgba('RoyalBlue')
mis_colores = matplotlib.colors.ListedColormap([gold, blue, purple, green],
'mis_colores')
In [246]:
plt.scatter(z_approx[:, 0], z_approx[:, 1], c=labels, s=40,
cmap=mis_colores, alpha=0.5, vmin = 0, vmax=3)
plt.scatter(centers[:, 0], centers[:, 1], c='red', s=200, alpha=0.5);
plt.xlabel('Componente 1 \n (intensificación)')
plt.ylabel('Componente 2 \n (área agríola)')
plt.title('Componentes calculadas con el método de Newton')
plt.show()
Comentario final sobre la nota:
Si bien el enfoque que se revisó en esta nota resaltó la alta dependencia del punto inicial para el método de Newton, representa un ejemplo de uso de métodos de descenso para resolver el problema de componentes principales. En Sparse Principal Component Analysis se utiliza un enfoque similar al de la función glmnet vía 4.2.Descenso_por_coordenadas_R. En este paper además de realizar la reducción de dimensionalidad, las componentes principales tienen una mejor interpretabilidad debido a que los loadings resultantes son cero y por tanto se seleccionan variables más importantes. Ver Sparse PCA para referencias de métodos que buscan resolver Sparse PCA.
Referencias:
S. P. Boyd, L. Vandenberghe, Convex Optimization, Cambridge University Press, 2009.