Por que PCA de dados por meio de SVD dos dados?

Esta pergunta é sobre uma maneira eficiente de calcular os componentes principais.

Muitos textos sobre PCA linear advogam o uso da decomposição de valor singular dos dados casewise . Ou seja, se temos dados $\bf X$ e queremos substituir as variáveis (suas colunas ) por componentes principais, fazemos SVD: $\bf X=USV'$ , valores singulares (raízes quadradas dos valores próprios) que ocupam a diagonal principal de $\bf S$ , os autovetores direitos $\bf V$ são a matriz de rotação ortogonal dos eixos-variáveis em eixos-componentes; os autovetores esquerdos $\bf U$ são como $\bf V$ , apenas para os casos. Podemos então calcular os valores dos componentes como $\bf C=XV=US$ .
Outra maneira de realizar o PCA de variáveis é via decomposição da matriz quadrada (isto é, pode ser correlações ou covariâncias etc., entre as variáveis). A decomposição pode ser uma decomposição autônoma ou decomposição de valor singular: com matriz semidefinida positiva simétrica quadrada, eles obterão o mesmo resultado com valores próprios na diagonal de e conforme descrito anteriormente. Os valores dos componentes serão . $\bf R=X'X$ $\bf R$ $\bf R=VLV'$ $\bf L$ $\bf V$ $\bf C=XV$

Agora, minha pergunta: se os dados são uma matriz grande e o número de casos é (geralmente um caso) muito maior que o número de variáveis, então o caminho (1) deve ser muito mais lento que o caminho (2), porque a maneira (1) aplica um algoritmo bastante caro (como SVD) a uma grande matriz; ele calcula e armazena uma enorme matriz que realmente não precisamos no nosso caso (o PCA de variáveis). Se sim, então por que tantos livros didáticos parecem advogar ou apenas mencionar a única maneira (1)? Talvez seja eficiente e estou perdendo alguma coisa? $\bf X$ $\bf U$

pca algorithms svd matrix-decomposition ttnphns
fonte

Geralmente, estamos interessados apenas nos poucos componentes principais que explicam a maior parte da variação. É possível fazer um SVD reduzido; por exemplo, se

é de dimensão

onde

seguida da função vai calcular apenas o primeiro

esquerdo e direito por vectores singulares padrão.

X

$X$

N \times p

$N \times p$

p << N

$p << N$ Rsvd

p

$p$

M. Berk

@ M.Berk:

, no entanto, é o mesmo em ambas as abordagens: elas produzem resultados equivalentes (iguais a alterações de sinal). Além disso, por exemplo, R calcula

somente se solicitado.

p

$p$

C

$\mathbf C$

cbeleites apoia Monica

Você tem uma referência para o caminho (1)? Só estou ciente de que o PCA está sendo implementado via SVD na matriz de covariância (ou seja, caminho 2), pois isso evita alguns problemas numéricos e se adapta obviamente à dimensionalidade, não ao tamanho do conjunto de dados. Maneira (1) Eu chamaria SVD, não PCA. Eu só vi isso em um contexto SVD puro, onde não seria possível realizar uma decomposição completa.

Anony-Mousse

@ Anony-Mousse, Apenas uma a mencionar, na Joliffe, Principal component analysis, 2nd ed.verdade, Joliffe descreve os dois lados, mas no capítulo principal do PCA ele diz apenas o modo 1, tanto quanto me lembro.

precisa saber é o seguinte

@ Anony-Mousse, o Caminho 1 para mim é importante do ponto de vista teórico, porque mostra claramente como o PCA está diretamente relacionado à simples análise de correspondência .

precisa saber é o seguinte

Respostas:

Aqui estão os meus 2ct no tópico

A aula de quimiometria onde eu aprendi o PCA usou a solução (2), mas não era orientada numericamente, e minha aula numérica era apenas uma introdução e não discutia SVD até onde me lembro.
Se eu entendo Holmes: SVD rápido para matrizes de grande escala corretamente, sua ideia foi usada para obter um SVD computacionalmente rápido de matrizes longas.
Isso significaria que uma boa implementação de SVD pode seguir internamente (2) se encontrar matrizes adequadas (não sei se ainda existem melhores possibilidades). Isso significaria que, para uma implementação de alto nível, é melhor usar o SVD (1) e deixar para o BLAS cuidar de qual algoritmo usar internamente.

Verificação prática rápida: o DVD do OpenBLAS parece não fazer essa distinção, em uma matriz de 5e4 x 100, svd (X, nu = 0)assume mediana de 3,5 s, enquanto svd (crossprod (X), nu = 0)leva 54 ms (chamado de R com microbenchmark).
A quadratura dos valores próprios, é claro, é rápida e, até o momento, os resultados das duas chamadas são 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

atualização: Dê uma olhada em Wu, W .; Massart, D. & de Jong, S .: Os algoritmos PCA do kernel para dados amplos. Parte I: Teoria e algoritmos, Chemometrics and Intelligent Laboratory Systems, 36, 165 - 172 (1997). DOI: http://dx.doi.org/10.1016/S0169-7439(97)00010-5

Este artigo discute propriedades numéricas e computacionais de 4 algoritmos diferentes para PCA: SVD, decomposição de eigen (EVD), NIPALS e POWER.

Eles estão relacionados da seguinte maneira:

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

O contexto do trabalho é amplo e funciona em (kernel PCA) - essa é exatamente a situação oposta à que você pergunta. Portanto, para responder sua pergunta sobre o comportamento da matriz longa, você precisa trocar o significado de "kernel" e "clássico". $\mathbf X^{(30 \times 500)}$ $\mathbf{XX'}$

performance comparison

Não é de surpreender que EVD e SVD mudem de lugar, dependendo se os algoritmos clássicos ou de kernel são usados. No contexto desta pergunta, isso significa que um ou outro pode ser melhor dependendo da forma da matriz.

Porém, a partir da discussão sobre SVD e EVD "clássicos", fica claro que a decomposição de é uma maneira muito usual de calcular o PCA. No entanto, eles não especificam qual algoritmo SVD é usado, exceto o uso da função 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

cbeleites suporta Monica
fonte

Portanto, seu teste (3,5 s vs 54 ms) suporta minha linha de que "o caminho 1" é consideravelmente mais lento. Direita?

precisa saber é o seguinte

@ttnphns: sim. Mas como o svd é fornecido pelo BLAS, pode ser diferente com um BLAS diferente. Eu esperava que um bom BLAS otimizado fizesse algo assim. No entanto, não parece ser o caso do OpenBLAS. Estou com preguiça de verificar outros BLASs, mas talvez algumas pessoas possam verificar seus outros BLASs, para descobrirmos quais são otimizados para este caso e quais não são. (I enviado o desenvolvedor OpenBLAS e enviou-lhe um link para esta pergunta, talvez por isso ele pode acrescentar outras informações, por exemplo, razões para não mudar o algoritmo para o svd (X'X)para matrizes de comprimento.)

cbeleites suportes Monica

X^{'}

$X'$

X

$X$

n < p

$n<p$

X^{'} X

$X'X$

u_{n + 1} = X^{'} X u_{n} / | | X^{'} X u_{n} | |

$u_{n+1} = X'Xu_n / ||X'Xu_n||$

v_{1}

$v_1$

X

$X$

X^{'} X

$X'X$

X^{'} \times (X u_{n})

$X' \times (X u_n)$

X X^{T}

$\mathbf X \mathbf X^T$

Eu estava falando sobre sua atualização, onde a Nipals está envolvida. Confirmo que o Nipals não está envolvido no SVD de Lapack. Sobre seu experimento de benchmark, algo como também microbenchmark(X <- matrix(rnorm(5e6), ncol=100), Y <- t(X), svd(X), svd(Y), control=list(order="inorder"), times = 5)pode ser interessante.

Elvis

O SVD é mais lento, mas geralmente é considerado o método preferido devido à sua maior precisão numérica.

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

Aqui está o que está escrito na pca()função de ajuda do MATLAB :

Algoritmo de componente principal pcausado para executar a análise de componente [...] principal:

'svd' - padrão. Decomposição do valor singular (SVD) de X.

$n$ $p$

A última frase destaca o compromisso crucial de precisão e velocidade que está em jogo aqui.

$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$ (quarta linha) será o caminho mais rápido. O SVD da própria matriz de dados será mais lento de qualquer maneira.

$\mathbf X$ consigo pode levar a uma perda de precisão numérica. Aqui está um exemplo, adaptado da resposta da @ JM para Por que SVD na $X$ é preferível à composição automática de $XX^⊤$ no PCA em Math.SE.

Considere uma matriz de dados

X = (\begin{matrix} 1 & 1 & 1 \\ ϵ & 0 0 & 0 0 \\ 0 0 & ϵ & 0 0 \\ 0 0 & 0 0 & ϵ \end{matrix}),

$\mathbf X = \begin{pmatrix}1&1&1\\\epsilon & 0 & 0\\ 0 & \epsilon & 0 \\ 0 & 0 & \epsilon\end{pmatrix},$ às vezes chamada matriz de Läuchli (e omitamos a centralização deste exemplo). Seus valores singulares ao quadrado são

3 + ϵ^{2}

$3+\epsilon^2$ ,

ϵ^{2}

$\epsilon^2$ e

ϵ^{2}

$\epsilon^2$ . Levando

ϵ = 10^{- 5}

$\epsilon = 10^{-5}$ , podemos usar SVD e EIG para calcular estes valores:

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')')])

obtenção de resultados idênticos:

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

Mas agora $\epsilon = 10^{-10}$ podemos observar como o SVD ainda funciona bem, mas o EIG se decompõe:

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

O que acontece aqui é que o próprio cálculo da matriz de covariância quadratura o número da condição de $\mathbf X$ , especialmente no caso em que $\mathbf X$ possui algumas colunas quase colineares (ou seja, alguns valores singulares muito pequenos), primeiro calculando a matriz de covariância e depois calculando sua composição automática resultará em perda de precisão em comparação com o SVD direto.

I should add that one is often happy to ignore this potential [tiny] loss of precision and rather use the faster method.

amoeba says Reinstate Monica
fonte

Also for the SVD you can write (or just use) the algorithm for

X^{T}

$\mathbf X^T$ which has the opposite runtime behaviour (because the decomposition has symmetry wrt. transposing

X

$\mathbf X$ ) and otherwise the same numeric behaviour. Numeric stability is an important point (+1) - though I guess it depends a lot the data whether that matters. E.g. I typically far too few cases and noisy measurements - so the stability of my models is typically not limited by the numeric stability of underlying SVDs but by the patients/cell culture batches or instrumental noise.

cbeleites supports Monica

Thanks for the answer and for thorough consideration of pros and cons.

ttnphns

amoeba, may that be that you find time to show a concrete example where numerical stability suffers by eig() approach? (Readers will benefit: there is a point of trade-off between speed and stability. How one may decide in a concrete practical situation?)

ttnphns

@ttnphns I rewrote the whole answer, providing a concrete example. Take a look.

amoeba says Reinstate Monica

@amoeba, thank you very much for coming back and giving an example! I tried both epsilon examples in SPSS and got results like your's except for the very last line: instead of 3 0 -3.3307e-16 eigen in spss returned me 3 0 0. It looks as if the function has some in-built and fixed tolerance value beyond which it zero-offs. In this example, the function appeared as if to hack the knot of numerical instability by zeroing both tiny eigenvalues, the "0" and the "-16".

ttnphns