Dlaczego PCA danych za pomocą SVD danych?

To pytanie dotyczy skutecznego sposobu obliczania głównych składników.

Wiele tekstów na temat liniowego PCA opowiada się za dekompozycją danych w liczbie pojedynczej . Oznacza to, że jeśli mamy dane i chcemy zastąpić zmienne (jego kolumny ) głównymi składnikami, wykonujemy SVD: , wartości osobliwe (pierwiastki kwadratowe wartości własnych) zajmujące główną przekątną , prawe wektory własne są ortogonalną macierzą obrotu zmiennych osi w komponenty osi, lewe wektory własne są podobne do , tylko dla przypadków. Następnie możemy obliczyć wartości składników jako . $\bf X$ $\bf X=USV'$ $\bf S$ $\bf V$ $\bf U$ $\bf V$ $\bf C=XV=US$
Innym sposobem wykonania PCA zmiennych jest dekompozycja macierzy kwadratowej $\bf R=X'X$ (tzn. $\bf R$ może być korelacją lub kowariancją itp. Między zmiennymi). Rozkład może być rozkładem własnym lub rozkładem liczby pojedynczej: przy kwadratowej symetrycznej dodatniej macierzy półfinałowej dadzą ten sam wynik $\bf R=VLV'$ z wartościami własnymi jak przekątna $\bf L$ , i $\bf V$ jak opisano wcześniej. Wartości składników będą wynosić $\bf C=XV$ .

Moje pytanie: jeśli data $\bf X$ jest dużą macierzą, a liczba przypadków jest (co często jest przypadkiem) znacznie większa niż liczba zmiennych, to oczekuje się, że sposób (1) będzie znacznie wolniejszy niż sposób (2 ), ponieważ sposób (1) stosuje dość drogi algorytm (taki jak SVD) do dużej matrycy; oblicza i przechowuje ogromną macierz $\bf U$ której tak naprawdę nie potrzebujemy w naszym przypadku (PCA zmiennych). Jeśli tak, to dlaczego tak wiele podręczników wydaje się popierać lub po prostu wymieniać tylko sposób (1)? Może jest wydajny i czegoś mi brakuje?

— ttnphns
źródło

Zasadniczo interesuje nas tylko kilka głównych składników, które wyjaśniają większość wariancji. Możliwe jest zmniejszenie SVD; Na przykład, jeżeli jest wymiaru , gdzie wtedy jest funkcja obliczyć tylko pierwsza lewej i prawej wektorów osobliwych domyślnie.

X

$X$

N \times p

$N \times p$

p << N

$p << N$ Rsvd

p

$p$

— M. Berk,

@ M.Berk:

jednak taki sam w obu podejściach: dają równoważne wyniki (równe zmianom do znaku). Również np. R oblicza

tylko, jeśli jest o to poproszony.

p

$p$

C

$\mathbf C$

— cbeleites obsługuje Monikę

Czy masz referencje dotyczące sposobu (1)? Wiem tylko, że PCA jest implementowane za pomocą SVD na macierzy kowariancji (tj. Sposób 2), ponieważ pozwala to uniknąć pewnych problemów numerycznych i oczywiście skaluje się z wymiarowością, a nie z zestawem danych. Sposób (1) Nazwałbym SVD, a nie PCA w ogóle. Widziałem to tylko w czystym kontekście SVD, gdzie nie dokonałby się kompletnego rozkładu.

— Anony-Mousse

@ Anony-Mousse, Wystarczy wspomnieć, w Joliffe, Principal component analysis, 2nd ed.rzeczywistości Joliffe opisuje oba sposoby, ale w głównym rozdziale na temat PCA mówi o sposobie 1, o ile pamiętam.

— ttnphns

@ Anony-Mousse, Way 1 jest dla mnie ważny z teoretycznego punktu widzenia, ponieważ wyraźnie pokazuje, w jaki sposób PCA jest bezpośrednio związane z prostą analizą korespondencji .

— ttnphns

Odpowiedzi:

Oto moje 2ct na ten temat

Wykład chemometryczny, w którym po raz pierwszy nauczyłem się PCA, wykorzystywał rozwiązanie (2), ale nie był zorientowany numerycznie, a mój wykład liczbowy był jedynie wstępem i nie omawiałem SVD, o ile pamiętam.
Jeśli dobrze rozumiem Holmes: Fast SVD dla macierzy wielkoskalowych , twój pomysł został wykorzystany do uzyskania obliczeniowego szybkiego SVD długich matryc.
Oznaczałoby to, że dobra implementacja SVD może wewnętrznie następować (2), jeśli napotka odpowiednie matryce (nie wiem, czy są jeszcze lepsze możliwości). Oznaczałoby to, że w przypadku implementacji wysokiego poziomu lepiej jest użyć SVD (1) i pozostawić BLAS-owi sprawdzenie, jakiego algorytmu użyć wewnętrznie.

Szybka praktyczna kontrola: Wydaje się, że svd OpenBLAS nie robi tego rozróżnienia na matrycy 5e4 x 100, svd (X, nu = 0)przyjmuje medianę 3,5 s, a svd (crossprod (X), nu = 0)zajmuje 54 ms (wywoływane z R za pomocą microbenchmark).
Kwadrat wartości własnych jest oczywiście szybki, a do tego wyniki obu wywołań są równoważne.

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

aktualizacja: spójrz na Wu, W .; Massart, D. i de Jong, S .: Algorytmy jądra PCA dla szerokich danych. Część I: Teoria i algorytmy, Chemometrics and Intelligent Laboratory Systems, 36, 165 - 172 (1997). DOI: http://dx.doi.org/10.1016/S0169-7439(97)00010-5

W tym artykule omówiono liczbowe i obliczeniowe właściwości 4 różnych algorytmów PCA: SVD, rozkład własny (EVD), NIPALS i POWER.

Są one powiązane w następujący sposób:

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

Kontekst artykułu to szeroki i działają one na (jądro PCA) - jest to sytuacja odwrotna do tej, o którą pytasz. Aby odpowiedzieć na pytanie dotyczące długiego zachowania macierzy, musisz wymienić znaczenie „jądra” i „klasycznego”. $\mathbf X^{(30 \times 500)}$ $\mathbf{XX'}$

porównanie wydajności

Nic dziwnego, że EVD i SVD zmieniają miejsca w zależności od tego, czy stosowane są algorytmy klasyczne, czy jądra. W kontekście tego pytania oznacza to, że jedno lub drugie może być lepsze w zależności od kształtu matrycy.

Ale z ich dyskusji na temat „klasycznego” SVD i EVD jasno wynika, że rozkład jest bardzo powszechnym sposobem obliczania PCA. Nie określają jednak, który algorytm SVD jest używany, poza tym, że używają funkcji Matlaba . $\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 obsługuje Monikę
źródło

Tak więc twoje testy (3,5 s vs 54 msek) obsługują moją linię, że „droga 1” jest znacznie wolniejsza. Dobrze?

— ttnphns,

@ttnphns: tak. Ale jak svd jest dostarczany przez BLAS, który może być inny z innym BLAS. Spodziewałbym się, że dobry zoptymalizowany BLAS robi coś takiego. Jednak wydaje się, że tak nie jest w przypadku OpenBLAS. Jestem zbyt leniwy, by sprawdzać inne BLASY, ale może kilka osób mogłoby sprawdzić inne BLASY, więc dowiemy się, które z nich są zoptymalizowane dla tego przypadku, a które nie. (Wysłałem e-mail do programisty OpenBLAS i wysłałem mu link do tego pytania, więc może może dodać jakieś informacje, np. Powody, dla których nie przestawił algorytmu svd (X'X)na długie matryce.)

— cbeleites obsługuje Monikę

Niektóre punkty wymagają wyjaśnienia (do mnie). Czy „metody jądra” można streścić jako „działające na

zamiast

gdy

”? jeśli tak, jest to dość trywialne. Nie znam MOCY, ale znam NIPALS, który oblicza wektory własne

przez iterację

(jego zbieżność z 1. wektorem własnym

, wtedy musisz zaktualizować

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$ obliczyć sekundę itp.). Istnieją dwa sposoby wykonania NIPALS, (1) możesz wstępnie obliczyć

lub (2) możesz wykonać produkt jako

, który sposób jest tutaj użyty? Wydaje mi się, że użyto (1), co może być niesprawiedliwe.

X^{'} X

$X'X$

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

$X' \times (X u_n)$

— Elvis

@Elvis: a) Metody jądra mają więcej niż tylko przetwarzanie na

, patrz np . Stats.stackexchange.com/questions/2499/.... Równoważność jest trywialna w przypadku PCA (nie ma znaczenia, czy zaczynasz od uzyskania wektorów prawych, czy lewych), ale nie w przypadku innych metod. b) „jaki sposób robić NIPALS” opiera się na tej samej ogólnej zasadzie. Wybór algorytmu SVD zależy od twojego BLAS, a właściwie sądzę, że NIPALS nie jest w to zaangażowany. Zauważ, że mój czas obejmuje obliczenie iloczynu krzyżowego.

X X^{T}

$\mathbf X \mathbf X^T$

— cbeleites wspiera Monikę

Mówiłem o twojej aktualizacji, w której uczestniczy Nipals. Potwierdzam, że Nipals nie bierze udziału w SVD Lapacka. Jeśli chodzi o eksperyment porównawczy, coś podobnego microbenchmark(X <- matrix(rnorm(5e6), ncol=100), Y <- t(X), svd(X), svd(Y), control=list(order="inorder"), times = 5)może być również interesujące.

— Elvis

SVD jest wolniejsze, ale często jest uważane za preferowaną metodę ze względu na wyższą dokładność numeryczną.

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

Oto, co napisano w pomocypca() funkcji MATLAB :

Algorytm głównego elementu, który pcawykorzystuje się do wykonania analizy głównego elementu [...]:

„svd” - domyślnie. Rozkład wartości osobliwych (SVD) dla X.

$n$ $p$

Ostatnie zdanie podkreśla zasadniczy kompromis między dokładnością prędkości, jaki ma tutaj miejsce.

$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^⊤$

X = (\begin{matrix} 1 & 1 & 1 \\ ϵ & 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},$

3 + ϵ^{2}

$3+\epsilon^2$

ϵ^{2}

$\epsilon^2$

ϵ^{2}

$\epsilon^2$

ϵ = 10^{- 5}

$\epsilon = 10^{-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')')])

uzyskiwanie identycznych wyników:

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

$\epsilon = 10^{-10}$

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

$\mathbf X$ $\mathbf X$

Powinienem dodać, że często z przyjemnością zignorujesz tę potencjalną [małą] utratę precyzji i raczej zastosuję szybszą metodę.

— ameba mówi Przywróć Monikę
źródło

X^{T}

$\mathbf X^T$

X

$\mathbf X$

Dziękujemy za odpowiedź i dokładne rozważenie zalet i wad.

— ttnphns

ameba, czy to może być czas na pokazanie konkretnego przykładu, w którym stabilność numeryczna cierpi z powodu eig()podejścia? (Czytelnicy skorzystają: istnieje punkt kompromisu między szybkością a stabilnością. Jak można zdecydować w konkretnej praktycznej sytuacji?)

— ttnphns

@ttnphns Przepisałem całą odpowiedź, podając konkretny przykład. Spójrz.

— ameba mówi Przywróć Monikę

@amoeba, dziękuję bardzo za powrót i podanie przykładu! Próbowałem obu przykładów epsilon w SPSS i uzyskałem wyniki takie jak twoje, z wyjątkiem ostatniej linii: zamiast 3 0 -3.3307e-16eigen w spss zwróciło mi 3 0 0. Wygląda na to, że funkcja ma jakąś wbudowaną i ustaloną wartość tolerancji, powyżej której zeruje. W tym przykładzie funkcja wyglądała tak, jakby zhakować węzeł niestabilności numerycznej poprzez wyzerowanie zarówno niewielkich wartości własnych, „0”, jak i „-16”.

— ttnphns