Studiuję PCA z kursu Andrew Ng Coursera i innych materiałów. W pierwszym zadaniu cs224n na kursie NLP w Stanford oraz w filmie wykładowym Andrew Ng dokonują dekompozycji wartości pojedynczej zamiast dekompozycji wektorów własnych macierzy kowariancji, a Ng twierdzi nawet, że SVD jest liczbowo bardziej stabilny niż skład eigend.

Z mojego zrozumienia, dla PCA powinniśmy wykonać SVD macierzy danych (m,n)wielkości, a nie macierzy kowariancji (n,n)wielkości. I rozkład własny wektora macierzy kowariancji.

Dlaczego robią SVD macierzy kowariancji, a nie macierzy danych?

ameba udzieliła już dobrej odpowiedzi w komentarzach, ale jeśli chcesz formalnego argumentu, oto on.

Rozkład macierzy liczbie pojedynczej to , gdzie kolumny są wektorami własnymi a przekątne wpisy są pierwiastkami kwadratowymi jego wartości własnych, tj. .A = U Σ V T V A T A Σ σ i i = $A$$A=U\Sigma V^T$$V$$A^TA$$\Sigma$$\sigma_{ii}=\sqrt{\lambda_i(A^TA)}$

Jak wiecie, głównymi składnikami są ortogonalne rzuty zmiennych na przestrzeń wektorów własnych empirycznej macierzy kowariancji . Wariancja składników jest podana przez jej wartości własne, .λi($\frac{1}{n-1}A^TA$$\lambda_i(\frac{1}{n-1}A^TA)$

Rozważ dowolną macierz kwadratową , i wektor taki, że . Następnieα ∈ R v B v = λ $B$$\alpha \in \mathbb R$$v$$Bv=\lambda v$

1. $B^kv=\lambda^kv$
2. $\lambda(\alpha B) = \alpha\lambda( B)$

Zdefiniujmy . SVD z obliczy składową elektroniczną aby uzyskać$S=\frac{1}{n-1}A^TA$$S$$S^TS=\frac{1}{(n-1)^2}A^TAA^TA$

1. wektory własne , które przez właściwość 1 są wektorami$(A^TA)^TA^TA=A^TAA^TA$$A^TA$
2. te pierwiastki kwadratowe o wartości własnych , która z właściwości 2, a następnie 1, a następnie 2 ponownie, to .$\frac{1}{(n-1)^2}A^TAA^TA$$\sqrt{\frac{1}{(n-1)^2} \lambda_i(A^TAA^TA)} = \sqrt{\frac{1}{(n-1)^2} \lambda_i^2(A^TA)} = \frac{1}{n-1}\lambda_i(A^TA) = \lambda_i(\frac{1}{n-1}A^TA)$

Voilà!

Jeśli chodzi o stabilność liczbową, należałoby dowiedzieć się, jakie są zastosowane alogrithmy. Jeśli jesteś gotów, sądzę, że są to procedury LAPACK używane przez numpy:

• numpy.linalg.eig
• numpy.linalg.svd

Aktualizacja: Jeśli chodzi o stabilność, wydaje się, że implementacja SVD wykorzystuje podejście dziel i zwyciężaj, podczas gdy w składzie eigend zastosowano prosty algorytm QR. Nie mogę uzyskać dostępu do niektórych istotnych dokumentów SIAM z mojej instytucji (cięcia w badaniach), ale znalazłem coś, co mogłoby poprzeć ocenę, że procedura SVD jest bardziej stabilna.

W

Nakatsukasa, Yuji i Nicholas J. Higham. „Stabilne i wydajne algorytmy podziału i zdobycia widma dla symetrycznego rozkładu wartości własnych i SVD.” SIAM Journal on Scientific Computing 35.3 (2013): A1325-A1349.

porównują stabilność różnych algorytmów wartości własnych i wydaje się, że podejście dziel i zwyciężaj (używają tego samego co numpy w jednym z eksperymentów!) jest bardziej stabilne niż algorytm QR. To, wraz z twierdzeniami gdzie indziej, że metody D&C są rzeczywiście bardziej stabilne, popiera wybór Ng.

@amoeba miał doskonałe odpowiedzi na pytania PCA, w tym na temat stosunku SVD do PCA. Odpowiadając na twoje dokładne pytanie, podniosę trzy punkty:

• matematycznie nie ma różnicy, czy PCA oblicza się bezpośrednio na macierzy danych, czy na jej macierzy kowariancji
• różnica wynika wyłącznie z precyzji liczbowej i złożoności. Zastosowanie zastosowania SVD bezpośrednio do macierzy danych jest liczbowo bardziej stabilne niż do macierzy kowariancji
• SVD można zastosować do macierzy kowariancji w celu wykonania PCA lub uzyskania wartości własnych, w rzeczywistości jest to moja ulubiona metoda rozwiązywania problemów własnych

Okazuje się, że SVD jest bardziej stabilny niż typowe procedury dekompozycji wartości własnej, szczególnie w przypadku uczenia maszynowego. W uczeniu maszynowym łatwo jest uzyskać bardzo kolinearne regresory. SVD działa lepiej w tych przypadkach.

Oto kod Pythona, aby pokazać punkt. Stworzyłem wysoce współliniową macierz danych, uzyskałem jej macierz kowariancji i próbowałem uzyskać wartości własne tej ostatniej. SVD nadal działa, podczas gdy zwykły rozkład własny nie udaje się w tym przypadku.

import numpy as np
import math
from numpy import linalg as LA

np.random.seed(1)

# create the highly collinear series
T = 1000
X = np.random.rand(T,2)
eps = 1e-11
X[:,1] = X[:,0] + eps*X[:,1]

C = np.cov(np.transpose(X))
print('Cov: ',C)

U, s, V = LA.svd(C)
print('SVDs: ',s)

w, v = LA.eig(C)
print('eigen vals: ',w)

Wydajność:

Cov:  [[ 0.08311516  0.08311516]
 [ 0.08311516  0.08311516]]
SVDs:  [  1.66230312e-01   5.66687522e-18]
eigen vals:  [ 0.          0.16623031]

Aktualizacja

Odpowiadając na komentarz Federico Poloni, oto kod z testami stabilności SVD vs Eig na 1000 losowych próbkach tej samej matrycy powyżej. W wielu przypadkach Eig wykazuje 0 małych wartości własnych, co prowadziłoby do osobliwości macierzy, a SVD nie robi tego tutaj. SVD jest około dwa razy bardziej precyzyjny przy niewielkim określaniu wartości własnej, co może, ale nie musi być ważne, w zależności od twojego problemu.

import numpy as np
import math
from scipy.linalg import toeplitz
from numpy import linalg as LA

np.random.seed(1)

# create the highly collinear series
T = 100
p = 2
eps = 1e-8

m = 1000 # simulations
err = np.ones((m,2)) # accuracy of small eig value
for j in range(m):
    u = np.random.rand(T,p)
    X = np.ones(u.shape)
    X[:,0] = u[:,0]
    for i in range(1,p):
        X[:,i] = eps*u[:,i]+u[:,0]

    C = np.cov(np.transpose(X))

    U, s, V = LA.svd(C)

    w, v = LA.eig(C)

    # true eigen values
    te = eps**2/2 * np.var(u[:,1])*(1-np.corrcoef(u,rowvar=False)[0,1]**2)
    err[j,0] = s[p-1] - te
    err[j,1] = np.amin(w) - te


print('Cov: ',C)
print('SVDs: ',s)
print('eigen vals: ',w)
print('true small eigenvals: ',te)

acc = np.mean(np.abs(err),axis=0)    
print("small eigenval, accuracy SVD, Eig: ",acc[0]/te,acc[1]/te)

Wydajność:

Cov:  [[ 0.09189421  0.09189421]
 [ 0.09189421  0.09189421]]
SVDs:  [ 0.18378843  0.        ]
eigen vals:  [  1.38777878e-17   1.83788428e-01]
true small eigenvals:  4.02633695086e-18
small eigenval, accuracy SVD, Eig:  2.43114702041 3.31970128319

Tutaj kod działa kod. Zamiast generować losową macierz kowariancji do testowania procedur, generuję macierz losowych danych z dwiema zmiennymi: gdzie - niezależne jednolite zmienne losowe. Zatem macierz kowariancji to gdzie - wariancje mundurów i współczynnik korelacji między im.u , v ( σ 2 1 σ 2

$$x_1=u\\ x_2=u+\varepsilon v$$ $u,v$ $$\begin{pmatrix} \sigma_1^2 & \sigma_1^2 + \varepsilon \rho \sigma_1 \sigma_2\\ \sigma_1^2 + \varepsilon \rho \sigma_1 \sigma_2 & \sigma_1^2 + 2 \varepsilon \rho \sigma_1 \sigma_2 + \varepsilon^2 \sigma_2^2\sigma^2\end{pmatrix}$$ $\sigma_1^2,\sigma_2^2,\rho$

Jego najmniejsza wartość własna: Mała wartość własna nie może być obliczona po prostu podłączając do formuły ze względu na ograniczoną precyzję, więc musisz go rozwinąć:

$$\lambda= \frac 1 2 \left(\sigma_2^2 \varepsilon^2 - \sqrt{\sigma_2^4 \varepsilon^4 + 4 \sigma_2^3 \rho \sigma_1 \varepsilon^3 + 8 \sigma_2^2 \rho^2 \sigma_1^2 \varepsilon^2 + 8 \sigma_2 \rho \sigma_1^3 \varepsilon + 4 \sigma_1^4} + 2 \sigma_2 \rho \sigma_1 \varepsilon + 2 \sigma_1^2\right)$$ $\varepsilon$ $$\lambda\approx \sigma_2^2 \varepsilon^2 (1-\rho^2)/2$$

Wykonuję symulacje realizacji macierzy danych, obliczam wartości własne symulowanej macierzy kowariancji i otrzymuję błędy .λ j e j = λ - λ$j=1,\dots,m$$\hat\lambda_j$$e_j=\lambda-\hat\lambda_j$

Użytkownikom Pythona chciałbym zauważyć, że w przypadku macierzy symetrycznych (takich jak macierz kowariancji) lepiej jest użyć numpy.linalg.eighfunkcji zamiast numpy.linalg.eigfunkcji ogólnej .

eighjest 9-10 razy szybszy niż eigna moim komputerze (niezależnie od rozmiaru matrycy) i ma lepszą dokładność (na podstawie testu dokładności @ Aksakal).

Nie jestem przekonany do wykazania korzyści z dokładności SVD przy małych wartościach własnych. Test Aksakala jest o 1-2 rzędy wielkości bardziej wrażliwy na losowy stan niż na algorytm (spróbuj wykreślić wszystkie błędy zamiast zmniejszać je do jednego absolutnego maksimum). Oznacza to, że małe błędy w macierzy kowariancji będą miały większy wpływ na dokładność niż wybór algorytmu składowego eigend. Nie ma to również związku z głównym pytaniem, które dotyczy PCA. Najmniejsze komponenty są ignorowane w PCA.

Podobny argument można postawić na temat stabilności liczbowej. Gdybym musiał użyć metody macierzy kowariancji dla PCA, rozłożyłbym ją za pomocą eighzamiast svd. Jeśli się nie powiedzie (czego tu jeszcze nie pokazano), prawdopodobnie warto przemyśleć problem, który próbujesz rozwiązać, zanim zaczniesz szukać lepszego algorytmu.

Aby odpowiedzieć na ostatnią część pytania: „Dlaczego robią SVD macierzy kowariancji, a nie macierzy danych?”. Uważam, że dzieje się tak ze względu na wydajność i pamięć. Zazwyczaj będzie bardzo dużą liczbą i nawet jeśli jest duże, spodziewalibyśmy się .n m ≫ $m$$n$$m \gg n$

Obliczenie macierzy kowariancji, a następnie wykonanie SVD na tym jest znacznie szybsze niż obliczenie SVD na pełnej macierzy danych w tych warunkach, dla tego samego wyniku.

Nawet przy dość małych wartościach wzrost wydajności jest współczynnikiem tysięcy (milisekund vs sekund). Uruchomiłem kilka testów na moim komputerze, aby porównać za pomocą Matlaba:

To tylko czas pracy procesora, ale potrzeby w zakresie pamięci masowej są tak samo ważne, jeśli nie większe. Jeśli spróbujesz SVD na matrycy milion na tysiąc w Matlabie, domyślnie wystąpi błąd, ponieważ potrzebuje działającej wielkości tablicy 7,4 TB.