Ttilizamos la base de datos Yale Face Database B recortada y reducida en resolución. El conjunto de datos comprende 2410 imágenes en escala de grises de 38 personas diferentes, recortadas y alineadas, y se puede cargar como:

utils::download.file('https://github.com/erichson/data/raw/master/R/faces.RData', 
                     destfile = 'faces.RData', 
                     method = 'auto')

Por conveniencia computacional, las imágenes de rostros de 96x84 se almacenan como vectores columna en la matriz de datos. Por ejemplo, el primer rostro se puede mostrar como:

load("faces.RData")
face <- matrix(rev(faces[ , 1]), nrow = 84, ncol = 96)
image(face, col = gray(0:255 / 255))

Para aproximar las k = 25 eigenfaces dominantes, podemos utilizar la función estándar de PCA en R:

faces.pca <- prcomp(t(faces), k = 25, center = TRUE, scale. = TRUE)
## Warning: In prcomp.default(t(faces), k = 25, center = TRUE, scale. = TRUE) :
##  extra argument 'k' will be disregarded

Nota que la matriz de datos debe ser transpuesta, de modo que cada columna corresponda a una ubicación de píxel en lugar de a una persona. Aquí, el análisis se realiza sobre la matriz de correlación estableciendo el argumento scale = TRUE. El rostro promedio se proporciona como el atributo center.

A continuación, utilizamos el paquete SPCA, pero establecemos los parámetros de ajuste en α = 0 y β = 0, lo que reduce el método a PCA:

library(sparsepca)
spca.results <- spca(t(faces), k=25, alpha=0, beta=0, center=TRUE, scale=TRUE)
## [1] "Iteration:    1, Objective: 9.71492e+05, Relative improvement Inf"
summary(spca.results) 
##                             PC1      PC2     PC3     PC4     PC5     PC6    PC7
## Explained variance     2901.541 2706.701 388.080 227.062 118.408 112.685 94.854
## Standard deviations      53.866   52.026  19.700  15.069  10.882  10.615  9.739
## Proportion of variance    0.360    0.336   0.048   0.028   0.015   0.014  0.012
## Cumulative proportion     0.360    0.695   0.744   0.772   0.786   0.800  0.812
##                           PC8    PC9   PC10   PC11   PC12   PC13   PC14   PC15
## Explained variance     84.035 66.930 64.939 55.089 48.498 44.253 40.652 37.802
## Standard deviations     9.167  8.181  8.058  7.422  6.964  6.652  6.376  6.148
## Proportion of variance  0.010  0.008  0.008  0.007  0.006  0.005  0.005  0.005
## Cumulative proportion   0.823  0.831  0.839  0.846  0.852  0.857  0.862  0.867
##                          PC16   PC17   PC18   PC19   PC20   PC21   PC22   PC23
## Explained variance     35.438 31.053 30.293 27.768 26.561 25.759 24.454 22.635
## Standard deviations     5.953  5.572  5.504  5.270  5.154  5.075  4.945  4.758
## Proportion of variance  0.004  0.004  0.004  0.003  0.003  0.003  0.003  0.003
## Cumulative proportion   0.871  0.875  0.879  0.882  0.886  0.889  0.892  0.895
##                          PC24   PC25
## Explained variance     21.771 20.188
## Standard deviations     4.666  4.493
## Proportion of variance  0.003  0.003
## Cumulative proportion   0.897  0.900

Las primeras 5 componentes principales explican aproximadamente el 79% de la variación total en los datos, mientras que las primeras 25 componentes principales explican más del 90%. Finalmente, los vectores propios pueden visualizarse como eigenfaces, por ejemplo, el primer vector propio (eigenface) se muestra de la siguiente manera:

layout(matrix(1:25, 5, 5, byrow = TRUE))
for(i in 1:25) {
  par(mar = c(0.5,0.5,0.5,0.5))
  img <- matrix(spca.results$loadings[,i], nrow=84, ncol=96)
  image(img[,96:1], col = gray((255:0)/255), axes=FALSE, xaxs="i", yaxs="i", xaxt='n', yaxt='n',ann=FALSE )
}

Las eigenfaces codifican tanto las características faciales holísticas como la iluminación.

Sin embargo, en muchas aplicaciones, es preferible obtener una representación sparse (dispersa). Esto se debe a que una representación dispersa es más fácil de interpretar. Además, ayuda a evitar el sobreajuste en entornos de datos de alta dimensión, donde el número de variables es mayor que el número de observaciones.

SPCA (Sparse Principal Component Analysis) intenta encontrar vectores de pesos dispersos (cargas), es decir, vectores propios aproximados con solo algunos valores activos (diferentes de cero). Podemos calcular los vectores propios dispersos de la siguiente manera:

USANDO EL GENERALIZED POWER METHOD

library(gpowerr)
# Transpose the data matrix
faces_t <- t(faces)

# Set the desired number of components
num_components <- 25

# Perform Sparse PCA
gpower_results <- gpower(faces_t, k = num_components, rho = 0.1)
summary(gpower_results) 
##                                   [,1]
## Proportion of Explained Variance 0.907
## Proportion of Sparsity           0.699
# Set up a 5x5 grid for visualization
layout(matrix(1:num_components, 5, 5, byrow = TRUE))

# Loop through each sparse principal component and plot it
for (i in 1:num_components) {
  par(mar = c(0.5, 0.5, 0.5, 0.5))
  
  # Reshape the sparse principal component vector into an image format
  img <- matrix(gpower_results$weights[, i], nrow = 84, ncol = 96)
  
  # Plot the grayscale image
  image(img[, 96:1], col = gray((255:0) / 255), axes = FALSE, xaxs = "i", yaxs = "i")
}

plot(1:length(gpower_results$comp_var), gpower_results$comp_var, 
     type = "o", pch = 16, col = "blue",
     xlab = "Principal Component", 
     ylab = "Variance Explained",
     main = "Variance Explained by Each Principal Component")

prop_var_explained <- gpower_results$comp_var / sum(gpower_results$comp_var) * 100

barplot(prop_var_explained, names.arg = 1:length(prop_var_explained),
        main = "Proportion of Variance Explained by Each Component",
        xlab = "Principal Component", ylab = "Percentage of Variance",
        col = "lightblue", border = "black", las = 2)

Ahora usando el robust sparse PCA

rspca.results <- rspca(t(faces), k=25, alpha=1e-4, beta=1e-1, verbose=1, max_iter=1000, tol=1e-4, center=TRUE, scale=TRUE)
## [1] "Iteration:    1, Objective: 8.59314e+06, Relative improvement Inf"
## [1] "Iteration:   11, Objective: 3.17759e+06, Relative improvement 4.97785e-02"
## [1] "Iteration:   21, Objective: 2.61205e+06, Relative improvement 6.56096e-03"
## [1] "Iteration:   31, Objective: 2.55024e+06, Relative improvement 7.58327e-04"
## [1] "Iteration:   41, Objective: 2.54254e+06, Relative improvement 1.24434e-04"

Los valores del objetivo para cada iteración se pueden graficar de la siguiente manera:

plot(log(rspca.results$objective), col='red', xlab='Number of iterations', ylab='Objective value')

El gráfico muestra la evolución del valor de la función objetivo a lo largo de las iteraciones del Sparse PCA (SPCA).

En este contexto, la función objetivo mide el error en la reconstrucción del conjunto de datos mediante componentes principales dispersos, con el objetivo de maximizar la varianza explicada mientras se impone una restricción de dispersión en los vectores propios. Un valor más alto de la función objetivo en las primeras iteraciones indica que la representación actual de los datos es ineficiente, ya que la varianza capturada es baja y el error de reconstrucción es alto.

A medida que avanzan las iteraciones, el valor de la función objetivo disminuye rápidamente, lo que sugiere que el modelo está encontrando una mejor representación de los datos con componentes más informativas y dispersas. Posteriormente, la reducción se desacelera y se estabiliza alrededor de las 30-40 iteraciones, indicando que el algoritmo ha convergido y ha encontrado una solución óptima, minimizando la función objetivo sin riesgo de sobreajuste

summary(rspca.results)
##                             PC1      PC2     PC3    PC4    PC5    PC6    PC7
## Explained variance     2616.248 2423.452 209.574 90.391 29.552 26.213 19.291
## Standard deviations      51.149   49.229  14.477  9.507  5.436  5.120  4.392
## Proportion of variance    0.324    0.301   0.026  0.011  0.004  0.003  0.002
## Cumulative proportion     0.324    0.625   0.651  0.662  0.666  0.669  0.671
##                           PC8   PC9  PC10  PC11  PC12  PC13  PC14  PC15  PC16
## Explained variance     14.769 9.607 8.984 6.472 4.814 3.855 3.085 2.670 2.343
## Standard deviations     3.843 3.100 2.997 2.544 2.194 1.963 1.756 1.634 1.531
## Proportion of variance  0.002 0.001 0.001 0.001 0.001 0.000 0.000 0.000 0.000
## Cumulative proportion   0.673 0.674 0.676 0.676 0.677 0.677 0.678 0.678 0.678
##                         PC17  PC18  PC19  PC20  PC21  PC22  PC23  PC24  PC25
## Explained variance     1.886 1.561 1.260 1.065 1.049 0.892 0.872 0.707 0.469
## Standard deviations    1.373 1.249 1.122 1.032 1.024 0.945 0.934 0.841 0.685
## Proportion of variance 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
## Cumulative proportion  0.679 0.679 0.679 0.679 0.679 0.679 0.680 0.680 0.680

Los primeros cinco componentes representan el 67% de la varianza

layout(matrix(1:25, 5, 5, byrow = TRUE))
for(i in 1:25) {
    par(mar = c(0.5,0.5,0.5,0.5))
    img <- matrix(rspca.results$loadings[,i], nrow=84, ncol=96)
    image(img[,96:1], col = gray((255:0)/255), axes=FALSE, xaxs="i", yaxs="i", xaxt='n', yaxt='n',ann=FALSE )
}

A diferencia de PCA, las cargas dispersas permiten contextualizar características localizadas. Si se desea, la solución puede hacerse aún más sparse aumentando el parámetro de ajuste α (alpha):

rspca.results <- rspca(t(faces), k=25, alpha=2e-4, beta=2e-1, verbose=1, max_iter=1000, tol=1e-4, center=TRUE, scale=TRUE)
## [1] "Iteration:    1, Objective: 1.43675e+07, Relative improvement Inf"
## [1] "Iteration:   11, Objective: 3.55885e+06, Relative improvement 2.66837e-02"
## [1] "Iteration:   21, Objective: 3.37726e+06, Relative improvement 4.72023e-04"
summary(rspca.results)
##                             PC1      PC2     PC3    PC4    PC5    PC6   PC7
## Explained variance     2378.526 2190.247 138.059 52.035 14.931 12.170 8.920
## Standard deviations      48.770   46.800  11.750  7.214  3.864  3.489 2.987
## Proportion of variance    0.295    0.272   0.017  0.006  0.002  0.002 0.001
## Cumulative proportion     0.295    0.567   0.584  0.590  0.592  0.593 0.595
##                          PC8   PC9  PC10  PC11  PC12  PC13  PC14  PC15  PC16
## Explained variance     6.333 3.913 3.382 2.479 1.574 1.181 0.846 0.705 0.626
## Standard deviations    2.517 1.978 1.839 1.574 1.255 1.087 0.920 0.840 0.791
## Proportion of variance 0.001 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
## Cumulative proportion  0.595 0.596 0.596 0.597 0.597 0.597 0.597 0.597 0.597
##                         PC17  PC18  PC19  PC20  PC21  PC22  PC23  PC24  PC25
## Explained variance     0.559 0.391 0.296 0.244 0.234 0.177 0.133 0.114 0.050
## Standard deviations    0.747 0.626 0.544 0.494 0.484 0.420 0.364 0.337 0.224
## Proportion of variance 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
## Cumulative proportion  0.597 0.597 0.597 0.597 0.597 0.597 0.597 0.597 0.597
layout(matrix(1:25, 5, 5, byrow = TRUE))
for(i in 1:25) {
    par(mar = c(0.5,0.5,0.5,0.5))
    img <- matrix(rspca.results$loadings[,i], nrow=84, ncol=96)
    image(img[,96:1], col = gray((255:0)/255), axes=FALSE, xaxs="i", yaxs="i", xaxt='n', yaxt='n',ann=FALSE )
}