Principal Components Analysis El objetivo principal de PCA es reducir la dimensionalidad. Que es la dimensionalidad y como lo hace?.
Sea la variable X^* una representacion de los datos en una menor dimension. El Objetivo es que si reconstruimos los datos al espacio original se pierda la menor cantidad de informacion.
Los datos se deben centrar antes de realizar cualquier projeccion.
$Y_{n,k} = (X_{n,d}-\mu_d) . W_{d,k}$
$X^*_{n,d} = Y_{n,k}.W^{-1}_{k,d}$
Como $W^{-1}_{k,d}$ es ortonormal, la inversa es la transpuesta $W^{-1}_{k,d}= W^T$.
$Min (X-X^*)^2$ , $W^T.W=I$
$=((X-\mu) - Y_{n,k}.W^{-1}_{k,d} )^2$
$=((X-\mu) - (X_{n,d}-\mu)W_{d,k}W^{-1}_{k,d})^2$
Sin embargo, minimizar el error, es similar a maximizar
$S = W^T \Sigma W$
$S=\sum_{i=1}^D{(X-\mu)(X-\mu)}^T$
$X_{n,d}- X_{n,k}^* = \sum_{i=1}^D{ W_i^T \Sigma W_i} - \sum_{i=1}^K{ W_i^T \Sigma W_i}$
Para reducir el error, las varianzas de los primeros k elementos van a ser mas grandes.
$det(A-\lambda)=0$
$m = [[1 2 4] [3 5 6] [1 0 2]]$
import numpy as np
import os
import pandas as pd
import matplotlib.pyplot as plt
import tarfile
import zipfile
from scipy.io import loadmat
from urllib.request import urlretrieve
from os.path import isfile, isdir
import seaborn as sns; sns.set()
%matplotlib inline
rng = np.random.RandomState(1)
X = np.dot(rng.rand(2, 2), rng.randn(2, 200)).T
plt.scatter(X[:, 0], X[:, 1])
plt.axis('equal');
m = np.mean(X,axis=0)
m
#calcular media
m = np.mean(X,axis=0)
X_center = X - m
X_center.mean(axis=0)
S = np.dot(X_center.T,X_center)
S.shape
S = np.dot(X_center.T,X_center)
S
# eigenvectors and eigenvalues for the from the scatter matrix
eig_val_sc, eig_vec_sc = np.linalg.eig(S)
eig_val_sc, eig_vec_sc
$W_{2,2}$
$W_{1,2} * X_{2,200}= X^*_{1,200}$
X_star = np.dot(X_center,eig_vec_sc[0,:].reshape(2,1))
Y = X_star *0
plt.scatter(X_star, Y, alpha=0.2)
from sklearn.decomposition import PCA
pca = PCA(n_components=2)
pca.fit(X)
def draw_vector(v0, v1, ax=None):
ax = ax or plt.gca()
arrowprops=dict(arrowstyle='->',
linewidth=2,
shrinkA=0,
shrinkB=0,
color='r')
ax.annotate('', v1, v0, arrowprops=arrowprops)
# plot data
plt.scatter(X[:, 0], X[:, 1], alpha=0.2)
for length, vector in zip(pca.explained_variance_, pca.components_):
v = vector * 3 * np.sqrt(length)
draw_vector(pca.mean_, pca.mean_ + v)
plt.axis('equal');
Vamos a utilizar las famosas Eigen faces, para mostrar como PCA puede ayudarnos en la practica
def download_files():
"""
Este metodo descarga los archivos de imagenes sino existen
"""
path_tar = os.path.join("data",'faces.zip')
if not isfile(path_tar):
urlretrieve(
'http://courses.media.mit.edu/2002fall/mas622j/proj/faces/rawdata.zip',
path_tar)
dest_path = "data/faces"
with zipfile.ZipFile(path_tar) as tar:
tar.extractall(dest_path)
tar.close()
download_files()
from scipy.io import loadmat
import glob
import random
import matplotlib.image as img
#Install pillow
from PIL import Image
d_name = 'data/faces/rawdata'
sample = 0.2
X = []
for nfile in os.listdir(d_name):
if random.random() <= sample:
bytes_read = open(os.path.join(d_name,nfile), "rb").read()
img = Image.frombytes('L', (128,128), bytes_read)
X.append(np.array(img).flatten())
X = np.array(X)
def plot_1_images(data, label="Image 1", ax=None):
fn_shape = lambda X: X.reshape(128,128)
fig = None
if ax is None:
fig, ax = plt.subplots(1,1, constrained_layout=True)
ax.imshow(fn_shape(data))
ax.set_title(label=label)
return fig,ax
def plot_3_images(data,ix_1, ix_2 , ix_3):
fn_shape = lambda X: X.reshape(128,128)
fig, ax = plt.subplots(1,3, constrained_layout=True)
ax[0].imshow(fn_shape(data[ix_1]))
ax[0].set_title(label="Image %s"% ix_1)
ax[1].imshow(fn_shape(data[ix_2]))
ax[1].set_title(label="Image %s"% ix_2)
ax[2].imshow(fn_shape(data[ix_2]))
ax[2].set_title(label="Image %s "% ix_3)
plt.show()
plot_3_images(X,1, 2 , 3)
X_mean = X.mean(axis=0)
X_center = X - X_mean
plot_3_images(X_center,1, 2 ,3)
plot_1_images(X_mean)
Calcular eigenvalues y eigen vectors
from sklearn.decomposition import PCA
pca = PCA(n_components=200)
# images x dim
X_reduced = pca.fit_transform(X_center)
eigen_values = pca.explained_variance_
eigen_faces = pca.components_
eigen_values
# Se puede ver como los primeras componentes contienen mas informacion
plt.plot(range(0,200), eigen_values/eigen_values[1])
plot_1_images(X[1])
Si pintamos las dimensiones, no tienen sentido. Porque cada elemento representa una mezcla de los pixeles de la foto original.
fig, ax = plt.subplots(1,1, constrained_layout=True)
ax.imshow(X_reduced[1].reshape(1,-1))
ax.set_title(label="Image 1")
plt.show()
print('Componentes', X_reduced.shape)
print('Eigenvectores', eigen_faces.shape)
face_1 = np.dot(X_reduced[1].reshape(1,-1),eigen_faces)
print(face_1.sum(axis=0))
eigen_faces[0].shape
face_c1 = np.dot(X_reduced[1].reshape(1,-1),eigen_faces)
face_c1.shape
plot_1_images(X_mean + face_1.sum(axis=0) )
Vemos que la reconstruccion de las imagenes funciono bien utilizando 200 eigen vectores.
face_10 = X_mean + np.sum(np.dot(X_reduced[1].reshape(1,-1)[:,:10],
eigen_faces[:10,:]),
axis=0)
face_50 = X_mean + np.sum(np.dot(X_reduced[1].reshape(1,-1)[:,:50],
eigen_faces[:50,:]),
axis=0)
face_100 = X_mean + np.sum(np.dot(X_reduced[1].reshape(1,-1)[:,:100],
eigen_faces[:100,:]),
axis=0)
face_200 = X_mean + np.sum(np.dot(X_reduced[1].reshape(1,-1)[:,:200],
eigen_faces[:200,:]),
axis=0)
fig, ax = plt.subplots(1,4,figsize=(15,15))
plot_1_images(face_20,ax=ax[0],label="Con 10 eigenfaces")
plot_1_images(face_20,ax=ax[1],label="Con 50 eigenfaces")
plot_1_images(face_100,ax=ax[2],label="Con 100 eigenfaces")
plot_1_images(face_200,ax=ax[3],label="Con 200 eigenfaces")
Como se puede ver en la imagen anterior. Las caracteristicas mas importantes son las formas de la cara. Despues las cejas y los ojos. Y por ultimo los detalles como el cabello, la boca y color de piel.
Despues de demostrar que 200 eigenvectors contienen la cantidad de informacion necesaria para reconstruir la imagen con calidad. Utilizaremos esas 200 dimensiones para encontrar las fotos mas similares.
min_dis = 100000000000000000
ix_1 = 0
ix_2 = 0
for i in range(200):
for j in range(i+1,200):
# excluir imagen en negro
if np.sum(X_reduced[i])>0:
dis_ij = sum(abs(X_reduced[i]-X_reduced[j]))
if dis_ij < min_dis:
min_dis = dis_ij
ix_1 = i
ix_2 = j
print("Imagenes similares son {0} y {1}".format(ix_1,ix_2))
plot_1_images(X[ix_1])
plot_1_images(X[ix_2])