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

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 $\bf X$ y queremos reemplazar las variables (sus columnas ) por componentes principales, hacemos SVD: $\bf 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 .$\bf S$$\bf V$$\bf U$$\bf V$$\bf C=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 .$\bf R=X'X$$\bf R$ $\bf R=VLV'$$\bf L$$\bf V$$\bf C=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?$\bf X$$\bf U$

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".$\mathbf X^{(30 \times 500)}$$\mathbf{XX'}$

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 .$\mathbf{X'X}$svd ()

    > 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

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

$\mathbf X$$\frac{1}{n-1}\mathbf X^\top \mathbf X$$\mathbf{XX}^\top$$n\ll p$

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.

$n$$p$

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

$1000\times 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.

$n \ll p$$\mathbf{XX}^\top$

$\mathbf X$$X$$XX^⊤$

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 $\epsilon = 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 de$\mathbf X$, especialmente en caso de que $\mathbf 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.