¿Por qué PCA de datos mediante SVD de los datos?

22

Esta pregunta trata sobre una manera eficiente de calcular los componentes principales.

  1. Muchos textos sobre PCA lineal abogan por el uso de la descomposición de valores singulares de los datos de casos . Es decir, si tenemos datos X y queremos reemplazar las variables (sus columnas ) por componentes principales, hacemos SVD: X=USV , valores singulares (raíces cuadradas de los valores propios) que ocupan la diagonal principal de , los vectores propios derechos son la matriz de rotación ortogonal de las variables de los ejes en componentes de los ejes, los vectores propios izquierdos son como , solo para los casos. Entonces podemos calcular los valores de los componentes como .SVUVC=XV=US

  2. Otra forma de hacer PCA de variables es mediante la descomposición de matriz cuadrada (es decir, pueden ser correlaciones o covarianzas , etc., entre las variables). La descomposición puede ser descomposición propia o descomposición de valores singulares: con una matriz semidefinida positiva simétrica cuadrada, darán el mismo resultado con valores propios que la diagonal de y como se describió anteriormente. Los valores de los componentes serán .R=XXR R=VLVLVdo=XV

Ahora, mi pregunta: si data es una matriz grande, y el número de casos es (que a menudo es un caso) mucho mayor que el número de variables, entonces se espera que way (1) sea ​​mucho más lento que way (2 ), porque way (1) aplica un algoritmo bastante costoso (como SVD) a una matriz grande; calcula y almacena una matriz enorme que realmente no necesitamos en nuestro caso (el PCA de las variables). Si es así, ¿ por qué tantos libros de texto parecen abogar o solo mencionar la única forma (1)? ¿Quizás es eficiente y me falta algo?XU

ttnphns
fuente
2
En general, solo nos interesan los pocos componentes principales que explican la mayor parte de la varianza. Es posible hacer una SVD reducida; por ejemplo, si es de dimensión N × p donde p < < N entonces 's función calcular sólo la primera p izquierda y vectores singulares derecha por defecto. XN×pp<<NRsvdp
M. Berk
1
@ M.Berk: , sin embargo, es el mismo en ambos enfoques: producen resultados equivalentes (iguales hasta los cambios de signo). Además, por ejemplo, R calcula C solo si se solicita. pC
cbeleites apoya a Monica el
¿Tiene referencias para way (1)? Solo soy consciente de que PCA se está implementando a través de SVD en la matriz de covarianza (es decir, la forma 2), ya que esto evita algunos problemas numéricos y obviamente se escala con la dimensionalidad, no con el tamaño del conjunto de datos. Manera (1) Llamaría SVD, no PCA en absoluto. Solo lo he visto en un contexto SVD puro, donde uno no realizaría una descomposición completa en realidad.
Anony-Mousse el
@ Anony-Mousse, solo uno para mencionar, en Joliffe, Principal component analysis, 2nd ed.realidad, Joliffe describe ambas formas, pero en el capítulo central sobre PCA dice sobre la manera 1, por lo que puedo recordar.
ttnphns
@ Anony-Mousse, el camino 1 para mí es importante desde el punto de vista teórico porque muestra claramente cómo PCA está directamente relacionado con el análisis de correspondencia simple .
ttnphns

Respuestas:

7

Aquí están mis 2ct sobre el tema

  • La conferencia de quimiometría donde aprendí por primera vez que PCA usaba la solución (2), pero no estaba orientada numéricamente, y mi conferencia de numérica fue solo una introducción y no discutí SVD hasta donde recuerdo.

  • Si entiendo Holmes: SVD rápida para matrices de gran escala correctamente, su idea se ha utilizado para obtener una SVD computacionalmente rápida de matrices largas.
    Eso significaría que una buena implementación de SVD puede seguir internamente (2) si encuentra matrices adecuadas (no sé si todavía hay mejores posibilidades). Esto significaría que para una implementación de alto nivel es mejor usar SVD (1) y dejar que BLAS se encargue de qué algoritmo usar internamente.

  • Comprobación práctica rápida: el svd de OpenBLAS no parece hacer esta distinción, en una matriz de 5e4 x 100, svd (X, nu = 0)toma una mediana de 3.5 s, mientras que svd (crossprod (X), nu = 0)toma 54 ms (llamado desde R con microbenchmark).
    La cuadratura de los valores propios, por supuesto, es rápida, y hasta eso los resultados de ambas llamadas son equivalentes.

    timing  <- microbenchmark (svd (X, nu = 0), svd (crossprod (X), nu = 0), times = 10)
    timing
    # Unit: milliseconds
    #                      expr        min         lq    median         uq        max neval
    #            svd(X, nu = 0) 3383.77710 3422.68455 3507.2597 3542.91083 3724.24130    10
    # svd(crossprod(X), nu = 0)   48.49297   50.16464   53.6881   56.28776   59.21218    10
    

actualización: Echa un vistazo a Wu, W .; Massart, D. y de Jong, S .: Los algoritmos de kernel PCA para datos amplios. Parte I: Teoría y algoritmos, Chemometrics and Intelligent Laboratory Systems, 36, 165 - 172 (1997). DOI: http://dx.doi.org/10.1016/S0169-7439(97)00010-5

Este artículo discute las propiedades numéricas y computacionales de 4 algoritmos diferentes para PCA: SVD, descomposición propia (EVD), NIPALS y POWER.

Están relacionados de la siguiente manera:

computes on      extract all PCs at once       sequential extraction    
X                SVD                           NIPALS    
X'X              EVD                           POWER

El contexto del documento es amplio , y funcionan en X X ' (kernel PCA): esta es la situación opuesta a la que usted pregunta. Entonces, para responder a su pregunta sobre el comportamiento de la matriz larga, debe intercambiar el significado de "núcleo" y "clásico".X(30×500)XX

comparación de rendimiento

No es sorprendente que EVD y SVD cambien de lugar dependiendo de si se usan los algoritmos clásico o del núcleo. En el contexto de esta pregunta, esto significa que uno u otro puede ser mejor dependiendo de la forma de la matriz.

Pero a partir de su discusión sobre SVD y EVD "clásica", está claro que la descomposición de es una forma muy habitual de calcular el PCA. Sin embargo, no especifican qué algoritmo SVD se usa más que el que usan la función de Matlab .XXsvd ()


    > sessionInfo ()
    R version 3.0.2 (2013-09-25)
    Platform: x86_64-pc-linux-gnu (64-bit)

    locale:
     [1] LC_CTYPE=de_DE.UTF-8       LC_NUMERIC=C               LC_TIME=de_DE.UTF-8        LC_COLLATE=de_DE.UTF-8     LC_MONETARY=de_DE.UTF-8   
     [6] LC_MESSAGES=de_DE.UTF-8    LC_PAPER=de_DE.UTF-8       LC_NAME=C                  LC_ADDRESS=C               LC_TELEPHONE=C            
    [11] LC_MEASUREMENT=de_DE.UTF-8 LC_IDENTIFICATION=C       

    attached base packages:
    [1] stats     graphics  grDevices utils     datasets  methods   base     

    other attached packages:
    [1] microbenchmark_1.3-0

loaded via a namespace (and not attached):
[1] tools_3.0.2

$ dpkg --list libopenblas*
[...]
ii  libopenblas-base              0.1alpha2.2-3                 Optimized BLAS (linear algebra) library based on GotoBLAS2
ii  libopenblas-dev               0.1alpha2.2-3                 Optimized BLAS (linear algebra) library based on GotoBLAS2
cbeleites apoya a Monica
fuente
Por lo tanto, su prueba (3.5 segundos frente a 54 ms) respalda mi línea de que "forma 1" es considerablemente más lenta. ¿Correcto?
ttnphns
1
@ttnphns: sí. Pero como el BLAS proporciona svd, eso podría ser diferente con un BLAS diferente. Esperaba que un buen BLAS optimizado haga algo como esto. Sin embargo, no parece ser el caso con OpenBLAS. Soy demasiado vago para verificar otros BLAS, pero tal vez algunas personas puedan verificar sus otros BLAS para descubrir cuáles están optimizados para este caso y cuáles no. (Le envié un correo electrónico al desarrollador de OpenBLAS y le envié un enlace a esta pregunta, para que tal vez pueda agregar información, por ejemplo, razones para no cambiar el algoritmo a svd (X'X)matrices largas).
cbeleites apoya a Monica el
XXnorte<pagsXXtunorte+1=XXtunorte/ /El |El |XXtunorteEl |El |v1XXXX×(Xtunorte)
XXT
Estaba hablando de tu actualización, donde Nipals está involucrado. Confirmo que Nipals no está involucrado en la SVD de Lapack. Acerca de su experimento de referencia, algo así también microbenchmark(X <- matrix(rnorm(5e6), ncol=100), Y <- t(X), svd(X), svd(Y), control=list(order="inorder"), times = 5)puede ser interesante.
Elvis
18

SVD es más lento, pero a menudo se considera el método preferido debido a su mayor precisión numérica.

X1norte-1XXXXnortepags

Esto es lo que está escrito en la ayuda de la pca()función de MATLAB :

Algoritmo del componente principal que se pcautiliza para realizar el análisis del componente principal [...]:

'svd': predeterminado. Descomposición del valor singular (SVD) de X.

nortepags

La última oración destaca la compensación crucial de velocidad y precisión que está en juego aquí.

1000×100

X = randn([1000 100]);

tic; svd(X); toc         %// Elapsed time is 0.004075 seconds.
tic; svd(X'); toc        %// Elapsed time is 0.011194 seconds.
tic; eig(X'*X); toc      %// Elapsed time is 0.001620 seconds.
tic; eig(X*X'); toc;     %// Elapsed time is 0.126723 seconds.

nortepagsXX

XXXX

X=(111ϵ0 00 00 0ϵ0 00 00 0ϵ),
3+ϵ2ϵ2ϵ2ϵ=10-5 5
eps = 1e-5;
X = [1 1 1; eye(3)*eps];
display(['Squared sing. values of X: ' num2str(sort(svd(X),'descend').^2')])
display(['Eigenvalues of X''*X:       ' num2str(sort(eig(X'*X),'descend')')])

obteniendo resultados idénticos:

Squared sing. values of X: 3       1e-10       1e-10
Eigenvalues of X'*X:       3       1e-10       1e-10

Pero tomando ahora ϵ=10-10 Podemos observar cómo SVD todavía funciona bien pero EIG se descompone:

Squared sing. values of X: 3       1e-20       1e-20
Eigenvalues of X'*X:       3           0 -3.3307e-16

Lo que sucede aquí es que el mismo cálculo de la matriz de covarianza cuadra el número de condición deX, especialmente en caso de que X tiene algunas columnas casi colineales (es decir, algunos valores singulares muy pequeños), primero calcula la matriz de covarianza y luego calcula su descomposición propia dará como resultado una pérdida de precisión en comparación con la SVD directa.

Debo agregar que a menudo uno se alegra de ignorar esta potencial [pequeña] pérdida de precisión y más bien usar el método más rápido.

ameba dice Reinstate Monica
fuente
1
También para el SVD puede escribir (o simplemente usar) el algoritmo para XT que tiene el comportamiento de tiempo de ejecución opuesto (porque la descomposición tiene simetría wrt. transposición X) y, de lo contrario, el mismo comportamiento numérico. La estabilidad numérica es un punto importante (+1), aunque supongo que depende mucho de los datos si eso importa. Por ejemplo, típicamente tengo muy pocos casos y mediciones ruidosas, por lo que la estabilidad de mis modelos generalmente no está limitada por la estabilidad numérica de las SVD subyacentes sino por los pacientes / lotes de cultivos celulares o el ruido instrumental.
cbeleites apoya a Monica
Gracias por la respuesta y por considerar detenidamente los pros y los contras.
ttnphns
ameba, ¿puede ser que encuentres tiempo para mostrar un ejemplo concreto donde la estabilidad numérica sufre por el eig()enfoque? (Los lectores se beneficiarán: hay un punto de compensación entre velocidad y estabilidad. ¿Cómo se puede decidir en una situación práctica concreta?)
ttnphns
@ttnphns Reescribí toda la respuesta, proporcionando un ejemplo concreto. Echar un vistazo.
ameba dice Reinstate Monica
1
@amoeba, ¡muchas gracias por volver y dar un ejemplo! Probé ambos ejemplos de epsilon en SPSS y obtuve resultados como los de usted, excepto la última línea: en lugar de 3 0 -3.3307e-16eigen en spss me devolvió 3 0 0. Parece que la función tiene un valor de tolerancia incorporado y fijo más allá del cual se pone a cero. En este ejemplo, la función parecía cortar el nudo de la inestabilidad numérica al poner a cero ambos pequeños valores propios, el "0" y el "-16".
ttnphns