Warum PCA von Daten mittels SVD der Daten?

In dieser Frage geht es um eine effiziente Methode zur Berechnung von Hauptkomponenten.

Viele Texte zur linearen PCA befürworten die Verwendung der Singulärwertzerlegung der fallweisen Daten . Das heißt, wenn wir Daten und wollen die Variablen (seine ersetzen Spalten ) von Hauptkomponenten, wir tun SVD: (. Sq Wurzeln des Eigenwerts), singuläre Werte die Hauptdiagonale von Besatzungs , rechte Eigenvektoren sind die orthogonale Rotationsmatrix von Achsenvariablen in Achsenkomponenten, linke Eigenvektoren sind wie , nur für Fälle. Wir können dann Komponentenwerte berechnen als $\bf X$ $\bf X=USV'$ $\bf S$ $\bf V$ $\bf U$ $\bf V$ $\bf C=XV=US$ .
Eine andere Möglichkeit, eine PCA von Variablen durchzuführen, besteht in der Zerlegung der quadratischen Matrix (dh kann Korrelationen oder Kovarianzen usw. zwischen den Variablen sein). Die Zerlegung kann eine Eigenzerlegung oder eine Singulärwertzerlegung sein: Bei einer quadratischen symmetrischen positiven semidefiniten Matrix ergeben sie dasselbe Ergebnis mit Eigenwerten als Diagonale von und wie zuvor beschrieben. Komponentenwerte werden . $\bf R=X'X$ $\bf R$ $\bf R=VLV'$ $\bf L$ $\bf V$ $\bf C=XV$

Nun meine Frage: Wenn Daten eine große Matrix sind und die Anzahl der Fälle (was häufig der Fall ist) viel größer als die Anzahl der Variablen ist, wird erwartet, dass Weg (1) viel langsamer ist als Weg (2), weil way (1) einen recht teuren Algorithmus (wie SVD) auf eine große Matrix anwendet; es berechnet und speichert eine riesige Matrix die wir in unserem Fall wirklich nicht brauchen (die PCA von Variablen). Wenn ja, warum sprechen sich dann so viele Lehrbücher dafür aus oder erwähnen nur so (1)? Vielleicht ist es ist effizient und ich bin etwas fehlt? $\bf X$ $\bf U$

pca algorithms svd matrix-decomposition ttnphns
quelle

Im Allgemeinen interessieren uns nur die wenigen Hauptkomponenten, die den größten Teil der Varianz erklären. Es ist möglich, eine reduzierte SVD durchzuführen. wenn zum Beispiel

ist die Dimension

, wo

dann ‚s Funktion nur die erste Berechnung wird

durch Standard linke und rechte Singulärvektoren.

X

$X$

N \times p

$N \times p$

p << N

$p << N$ Rsvd

p

$p$

M. Berk

@ M.Berk:

jedoch in beiden Ansätzen gleich: Sie liefern gleichwertige Ergebnisse (gleich bis zum Vorzeichenwechsel). Auch zB berechnet R

nur auf Anfrage.

p

$p$

C

$\mathbf C$

cbeleites unterstützt Monica

Haben Sie Referenzen für Weg (1)? Mir ist nur bekannt, dass PCA über SVD in der Kovarianzmatrix implementiert wird (dh Weg 2), da dies einige numerische Probleme vermeidet und offensichtlich mit der Dimensionalität skaliert, nicht mit der Datensatzgröße. Weg (1) Ich würde SVD anrufen, überhaupt nicht PCA. Ich habe es nur in einem reinen SVD-Kontext gesehen, in dem man eigentlich keine vollständige Zerlegung durchführen würde.

Anony-Mousse

@ Anony-Mousse, Nur eine zu erwähnen, Joliffe, Principal component analysis, 2nd ed.Eigentlich beschreibt Joliffe beide Möglichkeiten, aber im Kernkapitel über PCA sagt er, soweit ich mich erinnern kann, nur über Weg 1.

TTNPHNS

@ Anony-Mousse, Weg 1 ist für mich vom theoretischen Standpunkt aus wichtig, da er klar zeigt, wie PCA direkt mit der einfachen Korrespondenzanalyse zusammenhängt .

TTNPHNS

Antworten:

Hier sind meine 2ct zum Thema

Die Vorlesung über Chemometrie, in der ich PCA zum ersten Mal lernte, verwendete Lösung (2), war jedoch nicht numerisch orientiert, und meine Vorlesung über Numerik war nur eine Einführung und behandelte SVD nicht, soweit ich mich erinnere.
Wenn ich Holmes: Schnelle SVD für Großmatrizen richtig verstehe , wurde Ihre Idee verwendet, um eine rechnerisch schnelle SVD für Langmatrizen zu erhalten.
Das würde bedeuten, dass eine gute SVD-Implementierung intern folgen könnte (2), wenn sie auf geeignete Matrizen stößt (ich weiß nicht, ob es noch bessere Möglichkeiten gibt). Dies würde bedeuten, dass es für eine Implementierung auf hoher Ebene besser ist, die SVD (1) zu verwenden und sie dem BLAS zu überlassen, um zu prüfen, welcher Algorithmus intern verwendet werden soll.

Schneller Praxistest: OpenBLAS 'svd scheint diese Unterscheidung nicht zu treffen. Auf einer Matrix von 5e4 x 100 svd (X, nu = 0)dauert es durchschnittlich 3,5 s, während es svd (crossprod (X), nu = 0)54 ms dauert (von R mit aufgerufen microbenchmark).
Die Quadratur der Eigenwerte ist natürlich schnell, und bis dahin sind die Ergebnisse beider Aufrufe gleichwertig.

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

Update: Schauen Sie sich Wu, W. an; Massart, D. & de Jong, S .: Die Kernel-PCA-Algorithmen für breite Daten. Teil I: Theorie und Algorithmen, Chemometrie und intelligente Laborsysteme, 36, 165 - 172 (1997). DOI: http://dx.doi.org/10.1016/S0169-7439(97)00010-5

In diesem Artikel werden numerische und rechnerische Eigenschaften von 4 verschiedenen Algorithmen für PCA erörtert: SVD, Eigenzerlegung (EVD), NIPALS und POWER.

Sie hängen wie folgt zusammen:

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

Der Kontext des Papiers sind breit , und sie arbeiten auf(KernelPCA) - das ist genau das Gegenteil Situation wie die Sie fragen. Um also Ihre Frage zum Langmatrixverhalten zu beantworten, müssen Sie die Bedeutung von "Kernel" und "Klassik" austauschen. $\mathbf X^{(30 \times 500)}$ $\mathbf{XX'}$

Leistungsvergleich

Es ist nicht überraschend, dass EVD und SVD die Plätze wechseln, je nachdem, ob der klassische oder der Kernel-Algorithmus verwendet wird. Im Zusammenhang mit dieser Frage bedeutet dies, dass das eine oder andere je nach Form der Matrix besser sein kann.

Aber aus ihrer Diskussion über "klassische" SVD und EVD geht hervor, dass die Zersetzung von $\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 unterstützt Monica
quelle

Ihr Test (3,5 Sek. Gegenüber 54 Sek.) Unterstützt meine Aussage, dass "Weg 1" erheblich langsamer ist. Recht?

TTNPHNS

@ttnphns: ja. Aber da svd von der BLAS bereitgestellt wird, könnte das bei einer anderen BLAS anders sein. Ich hätte erwartet, dass ein gut optimierter BLAS so etwas macht. Dies scheint jedoch bei OpenBLAS nicht der Fall zu sein. Ich bin zu faul, um andere BLASs zu überprüfen, aber vielleicht könnten ein paar Leute ihre anderen BLASs überprüfen, damit wir herausfinden, welche für diesen Fall optimiert sind und welche nicht. (Ich mailte die OpenBLAS Entwickler und schickte ihm einen Link zu dieser Frage, so dass er vielleicht ein paar Informationen hinzufügen können, zB Gründe für die nicht den Algorithmus auf die Schalt svd (X'X)für lange Matrizen.)

cbeleites unterstützt 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$ zur Berechnung der Sekunde usw). Es gibt zwei Möglichkeiten, NIPALS durchzuführen: (1) Sie können vorberechnen

X^{'} X

$X'X$

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

$X' \times (X u_n)$

X X^{T}

$\mathbf X \mathbf X^T$

Ich habe über dein Update gesprochen, bei dem Nipals involviert ist. Ich bestätige, dass Nipals nicht an Lapacks SVD beteiligt ist. Über Ihr Benchmark-Experiment kann auch so etwas microbenchmark(X <- matrix(rnorm(5e6), ncol=100), Y <- t(X), svd(X), svd(Y), control=list(order="inorder"), times = 5)interessant sein.

Elvis

SVD ist langsamer, wird jedoch aufgrund der höheren numerischen Genauigkeit häufig als bevorzugte Methode angesehen.

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

Hier ist , was in dem von MATLAB geschrieben pca()Funktion Hilfe :

Hauptkomponentenalgorithmus, pcamit dem die Hauptkomponentenanalyse durchgeführt wird [...]:

'svd' - Standard. Singularwertzerlegung (SVD) von X.

$n$ $p$

Der letzte Satz unterstreicht den entscheidenden Kompromiss zwischen Geschwindigkeit und Genauigkeit, der hier zum Tragen kommt.

$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.

Der schnellste Weg führt in diesem Fall über die Kovarianzmatrix (dritte Zeile). Natürlich, wenn $n \ll p$ $\mathbf{XX}^\top$

$\mathbf X$ $X$ $XX^⊤$ in PCA on Math.SE zu ersetzen.

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}$ können wir mit SVD und EIG folgende Werte berechnen:

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

identische Ergebnisse erhalten:

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

Aber jetzt nehmen $\epsilon = 10^{-10}$ wenn wir , können wir beobachten, wie SVD immer noch gut abschneidet, aber EIG zusammenbricht:

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

Was hier passiert, ist, dass die Berechnung der Kovarianzmatrix die Bedingungszahl von quadriert $\mathbf X$ , also vor allem für den Fall, wenn $\mathbf X$ Bei einigen nahezu kollinearen Spalten (dh bei einigen sehr kleinen Singularwerten) führt die Berechnung der Kovarianzmatrix und die anschließende Berechnung der Zerlegung zu einem Genauigkeitsverlust im Vergleich zur direkten SVD.

Ich sollte hinzufügen, dass man diesen möglichen [winzigen] Präzisionsverlust oftmals gerne ignoriert und lieber die schnellere Methode anwendet.

Amöbe sagt Reinstate Monica
quelle

Auch für die SVD können Sie den Algorithmus für schreiben (oder nur verwenden)

X^{T}

$\mathbf X^T$ welches das entgegengesetzte Laufzeitverhalten hat (weil die Zerlegung Symmetrie bezüglich Transponierung hat)

X

$\mathbf X$ ) und ansonsten das gleiche numerische Verhalten. Die numerische Stabilität ist ein wichtiger Punkt (+1) - obwohl ich denke, dass es sehr von den Daten abhängt, ob dies wichtig ist. ZB habe ich in der Regel viel zu wenige Fälle und verrauschte Messungen - daher ist die Stabilität meiner Modelle in der Regel nicht durch die numerische Stabilität der zugrunde liegenden SVDs beschränkt, sondern durch die Patienten- / Zellkulturchargen oder das instrumentelle Rauschen.

cbeleites unterstützt Monica

Vielen Dank für die Antwort und für die gründliche Prüfung der Vor- und Nachteile.

TTNPHNS

Amöbe, kann das sein, dass Sie Zeit finden, ein konkretes Beispiel zu zeigen, bei dem die numerische Stabilität durch eig()Annäherung leidet ? (Die Leser werden davon profitieren: Es gibt einen Kompromiss zwischen Geschwindigkeit und Stabilität. Wie kann man in einer konkreten praktischen Situation entscheiden?)

ttnphns

@ttnphns Ich habe die gesamte Antwort mit einem konkreten Beispiel umgeschrieben. Schau mal.

Amöbe sagt Reinstate Monica

@amoeba, vielen Dank, dass du zurückgekommen bist und ein Beispiel gegeben hast! Ich habe beide epsilon-Beispiele in SPSS ausprobiert und mit Ausnahme der allerletzten Zeile ähnliche Ergebnisse erhalten: Anstelle von 3 0 -3.3307e-16eigen in spss wurde ich zurückgegeben 3 0 0. Es sieht so aus, als ob die Funktion einen eingebauten und festen Toleranzwert hat, ab dem sie nullt. In diesem Beispiel schien die Funktion den Knoten der numerischen Instabilität zu hacken, indem beide winzigen Eigenwerte, die "0" und die "-16", auf Null gesetzt wurden.

TTNPHNS