Bester PCA-Algorithmus für eine Vielzahl von Funktionen (> 10K)?

54

Ich habe dies zuvor bei StackOverflow gefragt, aber es scheint, als wäre es hier angemessener, da es auf SO keine Antworten gab. Es ist eine Art Schnittstelle zwischen Statistik und Programmierung.

Ich muss Code schreiben, um PCA (Principal Component Analysis) durchzuführen. Ich habe die bekannten Algorithmen durchgesehen und diese implementiert , was, soweit ich das beurteilen kann, dem NIPALS-Algorithmus entspricht. Es funktioniert gut, um die ersten 2-3 Hauptkomponenten zu finden, scheint dann aber sehr langsam zu konvergieren (in der Größenordnung von Hunderten bis Tausenden von Iterationen). Hier sind die Details von dem, was ich brauche:

  1. Der Algorithmus muss effizient sein, wenn eine große Anzahl von Merkmalen (in der Größenordnung von 10.000 bis 20.000) und Stichprobengrößen in der Größenordnung von wenigen Hundert verarbeitet werden.

  2. Es muss ohne eine anständige lineare Algebra / Matrix-Bibliothek einigermaßen implementierbar sein, da die Zielsprache D ist, für die es noch keine gibt, und selbst wenn dies der Fall ist, würde ich es vorziehen, es nicht als Abhängigkeit zum betreffenden Projekt hinzuzufügen .

Nebenbei bemerkt, scheint R im selben Datensatz alle Hauptkomponenten sehr schnell zu finden, verwendet jedoch eine Singularwertzerlegung, die ich nicht selbst codieren möchte.

dsimcha
quelle
2
Es gibt viele öffentliche SVD-Algorithmen. Siehe en.wikipedia.org/wiki/… . Können Sie eine davon nicht verwenden oder anpassen? Außerdem ist R Open Source und unter einer GPL-Lizenz. Warum also nicht den Algorithmus ausleihen, wenn er die Aufgabe erfüllt?
Rob Hyndman
@Rob: Ich möchte es vermeiden, praktisch eine lineare Algebra-Bibliothek zu schreiben, und ich möchte auch das Copyleft der GPL vermeiden. Außerdem habe ich mir vorher einige Teile des R-Quellcodes angesehen, die im Allgemeinen nicht gut lesbar sind.
Dsimcha
4
Vermisse ich etwas? Sie haben> 10K Features, aber <1K Samples? Dies bedeutet, dass die letzten 9K-Komponenten beliebig sind. Möchten Sie alle 1K der ersten Komponenten?
Shabbychef
2
Auf jeden Fall müssen Sie SVD unbedingt implementieren. Dank der umfangreichen numerischen linearen Algebra-Forschung können Sie jetzt aus einer Vielzahl von Methoden auswählen, je nachdem, wie groß / klein, dünn / dicht Ihre Matrix ist oder ob Sie möchten nur die Singularwerte oder den vollständigen Satz von Singularwerten und linken / rechten Singularvektoren. Die Algorithmen sind meiner Meinung nach nicht besonders schwer zu verstehen.
JM ist kein Statistiker
Können Sie uns sagen, warum Sie PCA machen wollen?
Robin Girard

Antworten:

27

Ich habe die randomisierte SVD implementiert, wie in "Halko, N., Martinsson, PG, Shkolnisky, Y. & Tygert, M. (2010) angegeben. Ein Algorithmus für die Hauptkomponentenanalyse großer Datensätze. Arxiv preprint arXiv: 1007.5510, 0526. Abgerufen am 1. April 2011 von http://arxiv.org/abs/1007.5510 . " Wenn Sie eine abgeschnittene SVD erhalten möchten, arbeitet diese sehr viel schneller als die svd-Variationen in MATLAB. Sie können es hier bekommen:

function [U,S,V] = fsvd(A, k, i, usePowerMethod)
% FSVD Fast Singular Value Decomposition 
% 
%   [U,S,V] = FSVD(A,k,i,usePowerMethod) computes the truncated singular
%   value decomposition of the input matrix A upto rank k using i levels of
%   Krylov method as given in [1], p. 3.
% 
%   If usePowerMethod is given as true, then only exponent i is used (i.e.
%   as power method). See [2] p.9, Randomized PCA algorithm for details.
% 
%   [1] Halko, N., Martinsson, P. G., Shkolnisky, Y., & Tygert, M. (2010).
%   An algorithm for the principal component analysis of large data sets.
%   Arxiv preprint arXiv:1007.5510, 0526. Retrieved April 1, 2011, from
%   http://arxiv.org/abs/1007.5510. 
%   
%   [2] Halko, N., Martinsson, P. G., & Tropp, J. A. (2009). Finding
%   structure with randomness: Probabilistic algorithms for constructing
%   approximate matrix decompositions. Arxiv preprint arXiv:0909.4061.
%   Retrieved April 1, 2011, from http://arxiv.org/abs/0909.4061.
% 
%   See also SVD.
% 
%   Copyright 2011 Ismail Ari, http://ismailari.com.

    if nargin < 3
        i = 1;
    end

    % Take (conjugate) transpose if necessary. It makes H smaller thus
    % leading the computations to be faster
    if size(A,1) < size(A,2)
        A = A';
        isTransposed = true;
    else
        isTransposed = false;
    end

    n = size(A,2);
    l = k + 2;

    % Form a real n×l matrix G whose entries are iid Gaussian r.v.s of zero
    % mean and unit variance
    G = randn(n,l);


    if nargin >= 4 && usePowerMethod
        % Use only the given exponent
        H = A*G;
        for j = 2:i+1
            H = A * (A'*H);
        end
    else
        % Compute the m×l matrices H^{(0)}, ..., H^{(i)}
        % Note that this is done implicitly in each iteration below.
        H = cell(1,i+1);
        H{1} = A*G;
        for j = 2:i+1
            H{j} = A * (A'*H{j-1});
        end

        % Form the m×((i+1)l) matrix H
        H = cell2mat(H);
    end

    % Using the pivoted QR-decomposiion, form a real m×((i+1)l) matrix Q
    % whose columns are orthonormal, s.t. there exists a real
    % ((i+1)l)×((i+1)l) matrix R for which H = QR.  
    % XXX: Buradaki column pivoting ile yapılmayan hali.
    [Q,~] = qr(H,0);

    % Compute the n×((i+1)l) product matrix T = A^T Q
    T = A'*Q;

    % Form an SVD of T
    [Vt, St, W] = svd(T,'econ');

    % Compute the m×((i+1)l) product matrix
    Ut = Q*W;

    % Retrieve the leftmost m×k block U of Ut, the leftmost n×k block V of
    % Vt, and the leftmost uppermost k×k block S of St. The product U S V^T
    % then approxiamtes A. 

    if isTransposed
        V = Ut(:,1:k);
        U = Vt(:,1:k);     
    else
        U = Ut(:,1:k);
        V = Vt(:,1:k);
    end
    S = St(1:k,1:k);
end

Erstellen Sie zum Testen einfach ein Bild in demselben Ordner (Sie können die Matrix auch als große Matrix selbst erstellen).

% Example code for fast SVD.

clc, clear

%% TRY ME
k = 10; % # dims
i = 2;  % # power
COMPUTE_SVD0 = true; % Comment out if you do not want to spend time with builtin SVD.

% A is the m×n matrix we want to decompose
A = im2double(rgb2gray(imread('test_image.jpg')))';

%% DO NOT MODIFY
if COMPUTE_SVD0
    tic
    % Compute SVD of A directly
    [U0, S0, V0] = svd(A,'econ');
    A0 = U0(:,1:k) * S0(1:k,1:k) * V0(:,1:k)';
    toc
    display(['SVD Error: ' num2str(compute_error(A,A0))])
    clear U0 S0 V0
end

% FSVD without power method
tic
[U1, S1, V1] = fsvd(A, k, i);
toc
A1 = U1 * S1 * V1';
display(['FSVD HYBRID Error: ' num2str(compute_error(A,A1))])
clear U1 S1 V1

% FSVD with power method
tic
[U2, S2, V2] = fsvd(A, k, i, true);
toc
A2 = U2 * S2 * V2';
display(['FSVD POWER Error: ' num2str(compute_error(A,A2))])
clear U2 S2 V2

subplot(2,2,1), imshow(A'), title('A (orig)')
if COMPUTE_SVD0, subplot(2,2,2), imshow(A0'), title('A0 (svd)'), end
subplot(2,2,3), imshow(A1'), title('A1 (fsvd hybrid)')
subplot(2,2,4), imshow(A2'), title('A2 (fsvd power)')

Schnelle SVD

Wenn ich es auf meinem Desktop für ein Image der Größe 635 * 483 ausführe, erhalte ich

Elapsed time is 0.110510 seconds.
SVD Error: 0.19132
Elapsed time is 0.017286 seconds.
FSVD HYBRID Error: 0.19142
Elapsed time is 0.006496 seconds.
FSVD POWER Error: 0.19206

Wie Sie sehen, ist es bei niedrigen Werten von kmehr als 10-mal schneller als bei Verwendung von Matlab SVD. Übrigens benötigen Sie möglicherweise die folgende einfache Funktion für die Testfunktion:

function e = compute_error(A, B)
% COMPUTE_ERROR Compute relative error between two arrays

    e = norm(A(:)-B(:)) / norm(A(:));
end

Ich habe die PCA-Methode nicht hinzugefügt, da die Implementierung mit SVD unkompliziert ist. Sie können diesen Link überprüfen , um ihre Beziehung zu sehen.

Petrichor
quelle
12

Sie könnten versuchen, mit ein paar Optionen.

1- Bestrafte Matrixzersetzung . Sie wenden einige Strafbedingungen auf die u und v an, um eine gewisse Sparsamkeit zu erzielen. Schneller Algorithmus, der für Genomdaten verwendet wurde

Siehe Whitten Tibshirani. Sie haben auch ein R-Päckchen. "Eine bestrafte Matrixzerlegung mit Anwendungen für spärliche Hauptkomponenten und kanonische Korrelationsanalyse."

2- Randomisierte SVD . Da SVD ein Master-Algorithmus ist, kann eine sehr schnelle Approximation wünschenswert sein, insbesondere für die explorative Analyse. Mit randomisierter SVD können Sie PCA für große Datasets durchführen.

Siehe Martinsson, Rokhlin und Tygert "Ein randomisierter Algorithmus zur Zerlegung von Matrizen". Tygert hat Code für eine sehr schnelle Implementierung von PCA.

Unten sehen Sie eine einfache Implementierung der randomisierten SVD in R.

ransvd = function(A, k=10, p=5) {
  n = nrow(A)
  y = A %*% matrix(rnorm(n * (k+p)), nrow=n)
  q = qr.Q(qr(y))
  b = t(q) %*% A
  svd = svd(b)
  list(u=q %*% svd$u, d=svd$d, v=svd$v)
}
pslice
quelle
+1 für bestrafte Matrixzerlegung. Dieses Paket ist ziemlich erstaunlich. Ich sollte aber wahrscheinlich erwähnen, dass es "Witten" heißt, falls die Leute Probleme haben, das Zitat zu finden. Schließlich sagte das OP, sie wollten nichts in R geschrieben haben, aber im Grunde genommen wird jedes große SVD-Paket ein C-, C ++ - oder Fortran-Backend haben, um die Geschwindigkeit zu erhöhen.
David J. Harris
3

Ich würde vorschlagen, Kernel-PCA zu testen, dessen zeitliche / räumliche Komplexität von der Anzahl der Beispiele (N) und nicht von der Anzahl der Features (P) abhängt. Kernel-PCA arbeitet grundsätzlich mit der NxN-Kernelmatrix (Matrix der Ähnlichkeiten zwischen den Datenpunkten) und nicht mit der PxP-Kovarianzmatrix, die für große P schwer zu handhaben ist. Ein weiterer Vorteil von Kernel-PCA ist, dass es nichtlineare Projektionen lernen kann auch, wenn Sie es mit einem geeigneten Kernel verwenden. Siehe dieses Papier über Kernel-PCA .

Ebenholz1
quelle
2

Ich erinnere mich, dass es möglich ist, PCA durchzuführen, indem man die Eigenzerlegung von X ^ TX anstelle von XX ^ T berechnet und dann transformiert, um die PCs zu erhalten. Ich kann mich jedoch nicht an die Details erinnern, aber es ist in Jolliffes (ausgezeichnetem) Buch und ich werde es nachschlagen, wenn ich das nächste Mal bei der Arbeit bin. Ich würde die linearen Algebra-Routinen von zB Numerical Methods in C transliterieren, anstatt irgendeinen anderen Algorithmus zu verwenden.

Dikran Beuteltier
quelle
5
Meine Güte ... Die Konstruktion der Kovarianzmatrix ist für SVD niemals der beste Weg. Ich habe ein Beispiel gezeigt, warum die explizite Bildung der Kovarianzmatrix bei math.SE keine gute Idee ist: math.stackexchange.com/questions/3869/3871#3871 .
JM ist kein Statistiker
1

Es gibt auch die Bootstrap-Methode von Fisher et al. , Die für mehrere hundert Proben mit hohen Abmessungen ausgelegt ist.

Die Hauptidee der Methode lautet: "Resampling ist eine Transformation mit geringer Dimension". Wenn Sie also nur wenige (mehrere Hundert) hochdimensionale Proben haben, können Sie nicht mehr Hauptkomponenten als die Anzahl Ihrer Proben erhalten. Es ist daher sinnvoll, die Stichproben als sparsame Grundlage zu betrachten, die Daten auf den linearen Unterraum zu projizieren, der von diesen Vektoren überspannt wird, und die PCA innerhalb dieses kleineren Unterraums zu berechnen. Sie enthalten auch weitere Informationen zum Umgang mit dem Fall, dass möglicherweise nicht alle Proben im Speicher abgelegt sind.

Kolya Ivankov
quelle
0

Siehe Sam Roweis 'Artikel, EM-Algorithmen für PCA und SPCA .

ars
quelle
Der Wikipedia-Algorithmus zitiert dies und entspricht diesem für den Fall, dass jeweils eine Hauptkomponente gefunden wird.
Dsimcha
OK, ich sehe den Link jetzt. Dies ist ein ziemlich einfacher Ansatz, und genau wie bei Wikipedia gibt es Fortschritte bei dieser Grundidee. Überlegungen zufolge müssen Sie sich jedoch mit Kompromissen auseinandersetzen (Konvergenz in diesem Fall). Ich frage mich, ob Sie hier die richtige Frage stellen. Gibt es wirklich keine guten Bindungen zu Linalg-Bibliotheken für D?
ars