¿Cómo realizar el ACP para datos de muy alta dimensionalidad?

Question

Para realizar el análisis de componentes principales (ACP), hay que restar las medias de cada columna de los datos, calcular la matriz de coeficientes de correlación y, a continuación, hallar los vectores propios y los valores propios. Bueno, más bien, esto es lo que hice para implementarlo en Python, salvo que sólo funciona con matrices pequeñas porque el método para hallar la matriz del coeficiente de correlación (corrcoef) no me permite usar un array con alta dimensionalidad. Como tengo que usarlo para imágenes, mi implementación actual no me ayuda mucho.

He leído que es posible simplemente tomar su matriz de datos $D$ y calcula $DD^\top/n$ en lugar de $D^\top D/n$ pero a mí no me funciona. Bueno, no estoy exactamente seguro de entender lo que significa, además del hecho de que se supone que es un $n \times n$ en lugar de $p\times p$ (en mi caso $p\gg n$ ). Leí sobre ello en los tutoriales de eigenfaces, pero ninguno de ellos parecía explicarlo de forma que realmente pudiera entenderlo.

En resumen, ¿existe una descripción algorítmica sencilla de este método para que pueda seguirlo?

Answer 1

Lo que estás haciendo ahora está cerca, pero necesitas asegurarte de que multiplicas los vectores propios de (data . data.T) / lines a la izquierda por data.T para obtener los vectores propios de (data.T . data) / lines . Esto se denomina a veces el "truco de la transposición".

Aquí tienes más detalles. Supongamos que tenemos una matriz $A$ en el que desea realizar el ACP; para simplificar, supongamos que las columnas de $A$ ya se han normalizado para que tengan media cero, por lo que sólo tenemos que calcular los vectores propios de la matriz de covarianza $A^T A$ .

Ahora bien $A$ es un $m \times n$ matriz, con $n >> m$ entonces $A^T A$ es muy grande $n \times n$ matriz. Así, en lugar de calcular los vectores propios de $A^T A$ nos gustaría calcular los vectores propios de la mucho más pequeña $m \times m$ matriz $A A^T$ -- suponiendo que podamos averiguar la relación entre ambos. Entonces, ¿cómo son los vectores propios de $A^T A$ relacionados con los vectores propios de $A A^T$ ?

Sea $v$ sea un vector propio de $A A^T$ con valor propio $\lambda$ . Entonces

• $AA^T v = \lambda v$
• $A^T(A A^T v) = A^T(\lambda v)$
• $(A^T A)(A^T v) = \lambda (A^T v)$

En otras palabras, si $v$ es un vector propio de $A A^T$ entonces $A^T v$ es un vector propio de $A^T A$ con el mismo valor propio. Por lo tanto, al realizar un ACP en $A$ en lugar de encontrar directamente los vectores propios de $A^T A$ (que puede ser muy costoso), es más fácil encontrar los vectores propios $v$ de $AA^T$ y luego multiplicarlos a la izquierda por $A^T$ para obtener los vectores propios $A^T v$ de $A^T A$ .

Answer 2

La forma más sencilla de realizar un ACP estándar es centrar las columnas de la matriz de datos (suponiendo que las columnas correspondan a variables diferentes) restando las medias de las columnas y, a continuación, realizar una SVD. Los vectores singulares izquierdos, multiplicados por el valor singular correspondiente, corresponden a los componentes principales (estimados). Los vectores singulares de la derecha corresponden a las direcciones de los componentes principales (estimados), que son los mismos que los vectores propios del ACP. Los valores singulares corresponden a las desviaciones estándar de los componentes principales (multiplicados por un factor de raíz n, donde n es el número de filas de su matriz de datos) - lo mismo que la raíz cuadrada de los valores propios dados por PCA.

Si desea realizar el ACP en la matriz de correlaciones, deberá normalizar las columnas de la matriz de datos antes de aplicar la SVD. Esto equivale a restar las medias (centrado) y luego dividir por las desviaciones estándar (escalado).

Este será el enfoque más eficaz si desea el PCA completo. Puede verificar con un poco de álgebra que esto le da la misma respuesta que hacer la descomposición espectral de la matriz de covarianza de la muestra.

También existen métodos eficaces para calcular una SVD parcial, cuando sólo se necesitan algunos de los PC. Algunos de ellos son variantes de la iteración de potencia. En Algoritmo de Lanczos es un ejemplo que también está relacionado con los mínimos cuadrados parciales. Si tu matriz es enorme, puede que te vaya mejor un método aproximado. También hay razones estadísticas para regularizar el ACP cuando éste es el caso.

Answer 3

Parece que lo que quieres es el algoritmo NIPALS para realizar PCA. Es un algoritmo muy popular entre los estadísticos. Tiene muchas ventajas:

• Computacionalmente menos costoso que el SVD o los métodos de descomposición de valores propios si sólo se necesitan los primeros componentes.
• Tiene requisitos de almacenamiento más modestos en general porque la matriz de covarianza nunca se forma. Se trata de una propiedad muy importante para conjuntos de datos muy grandes.
• Puede manejar los datos que faltan en el conjunto de datos (aunque eso no es un problema en su problema, ya que se trata de imágenes).

Descripción
http://en.wikipedia.org/wiki/Non-linear_iterative_partial_least_squares

Algoritmo
He aquí una descripción sencilla y excelente del algoritmo (en la sección 1.2)
http://stats4.eng.mcmaster.ca/w/mediafiles/mediawiki/f/f7/Section-Extra-Class-1.pdf

Recuerde centrar primero la escala media antes de realizar el PCA, ya que es sensible a la escala.

Answer 4

Para añadir algo a la respuesta de Gilead, se trata de algoritmos computacionalmente menos costosos para los ACP truncados. NIPALS es de hecho muy popular, pero he tenido mucho éxito con métodos aproximados que realizan una sucesión de ajustes en datos parciales (lo que a menudo se llama PCA por proyección aleatoria). Esto se discutió en un metaoptimizar hilo.

Ya que mencionas Python, permíteme señalar que el algoritmo se implementa en el programa scikit-learn : el PCA clase. En concreto, se utiliza en un ejemplo que demuestra eigenfaces .

Answer 5

Si a alguien le sirve de ayuda, además del post de @raegtin, el siguiente código de R demuestra cómo se puede hacer PCA de alta dimensión de forma más rápida y eficiente.

Para la siguiente matriz $A_{m\times n}$ con $n>>m$ la matriz $A^TA$ es un $n\times n$ y no una matriz de rango completo, con rango $=m$ para que sólo el primer $m$ son distintos de cero.

Mientras que $AA^T$ es un $m \times m$ que es una matriz de rango completo, de nuevo con rango $m$ y $A^Tv$ es un vector propio de $A^TA$ cuando $v$ es un vector propio de $AA^T$ . Pero hay que tener en cuenta que $A^Tv$ calculados de esta manera no serán de longitud unitaria, por lo que tenemos que normalizar los vectores propios calculados de esta manera.

m <- 10
n <- 1000
A <- matrix(rnorm(m*n), m, n)

AtA <- t(A) %*% A
dim(AtA)
# [1] 1000 1000
qr(AtA)$rank
# [1] 10

AAt <- A %*% t(A)
dim(AAt)
# [1] 10 10
qr(AAt)$rank
# [1] 10

# eigenvalues and eigenvectors
res <- eigen(AtA)
val1 <- res$values
vec1 <- res$vectors

res <- eigen(AAt)
val2 <- res$values
vec2 <- res$vectors
vec2 <- t(A) %*% vec2

# eigenvalues
length(val1)
# [1] 1000
length(val2)
# [1] 10
val1 <- round(val1, 3)
val2 <- round(val2, 3)
val1
# [1] 1154.743 1098.774 1075.320 1049.960 1035.392  996.713  943.589  924.477 886.510  814.518 
# 0.000    0.000    0.000    0.000    0.000   0.000 ...
val2
# [1] 1154.743 1098.774 1075.320 1049.960 1035.392  996.713  943.589  924.477  886.510  814.518

all.equal(val1[1:10], val2)
# [1] TRUE   

# eigenvectors
sqrt(sum(vec1[,1]^2))
# [1] 1    
sqrt(sum(vec2[,1]^2))
# [1] 33.9815

# normalize eigenvectors
vec2 <- apply(vec2, 2, function(x) x / sqrt(sum(x^2)))
sqrt(sum(vec2[,1]^2))
# [1] 1

# compare run times
library(microbenchmark)
library(ggplot2)

mbm <- microbenchmark("PCA_AtA" = { 
    res <- eigen(AtA)
    val1 <- res$values
    vec1 <- res$vectors
  },
  "PCA_AAt" = {
    res <- eigen(AAt)
    val2 <- res$values
    vec2 <- res$vectors
    vec2 <- t(A) %*% vec2
  })

autoplot(mbm)