Najlepszy algorytm PCA dla ogromnej liczby funkcji (> 10 KB)?

54

Wcześniej zapytałem o to na StackOverflow, ale wydaje się, że może być bardziej odpowiednie tutaj, biorąc pod uwagę, że nie otrzymało żadnych odpowiedzi na SO. To trochę na styku statystyki i programowania.

Muszę napisać kod, aby wykonać PCA (Principal Component Analysis). Przejrzałem dobrze znane algorytmy i zaimplementowałem ten , który, o ile wiem, jest równoważny algorytmowi NIPALS. Działa dobrze do znajdowania pierwszych 2-3 głównych składników, ale wydaje się, że ich konwergencja staje się bardzo powolna (rzędu setek do tysięcy iteracji). Oto szczegóły tego, czego potrzebuję:

  1. Algorytm musi być wydajny w przypadku ogromnej liczby funkcji (od 10 000 do 20 000) i wielkości próbek rzędu kilkuset.

  2. Musi być racjonalnie możliwy do wdrożenia bez przyzwoitej biblioteki algebry liniowej / macierzy, ponieważ językiem docelowym jest D, który jeszcze go nie ma, a nawet gdyby tak było, wolałbym nie dodawać go jako zależności do danego projektu .

Na marginesie, wydaje się, że w tym samym zbiorze danych R wydaje się, że wszystkie główne składniki są bardzo szybkie, ale wykorzystuje rozkład pojedynczej wartości, co nie jest czymś, co chcę sam kodować.

dsimcha
źródło
2
Istnieje wiele publicznych algorytmów SVD. Zobacz en.wikipedia.org/wiki/… . Nie możesz użyć lub dostosować jednego z nich? Ponadto R jest oprogramowaniem typu open source i na licencji GPL, więc dlaczego nie pożyczyć jego algorytmu, jeśli spełnia swoje zadanie?
Rob Hyndman,
@Rob: Chciałbym uniknąć praktycznie pisania biblioteki algebry liniowej, a także chcę uniknąć copyleft GPL. Ponadto patrzyłem wcześniej na fragmenty kodu źródłowego R i ogólnie nie jest on zbyt czytelny.
dsimcha
4
Czy coś brakuje? Masz> 10 000 funkcji, ale <1 000 próbek? Oznacza to, że ostatnie komponenty 9K są dowolne. Czy chcesz mieć wszystkie 1K pierwszych elementów?
shabbychef
2
W każdym razie nie można uniknąć konieczności implementacji SVD, choć dzięki licznym badaniom z zakresu algebry liniowej istnieje teraz wiele metod do wyboru, w zależności od tego, jak duża / mała, rzadka / gęsta jest twoja matryca, lub jeśli chcesz tylko wartości w liczbie pojedynczej lub pełnego zestawu wartości w liczbie pojedynczej i wektorów w lewo / prawo w liczbie pojedynczej. Algorytmy nie są strasznie trudne do zrozumienia IMHO.
JM nie jest statystykiem
Czy możesz nam powiedzieć, dlaczego chcesz robić PCA?
robin girard

Odpowiedzi:

27

Zaimplementowałem Randomized SVD, jak podano w „Halko, N., Martinsson, PG, Shkolnisky, Y., i Tygert, M. (2010). Algorytm analizy głównej składowej dużych zbiorów danych. Arxiv preprint arXiv: 1007.5510, 0526. Źródło: 01 kwietnia 2011, od http://arxiv.org/abs/1007.5510 . ". Jeśli chcesz obciąć SVD, to naprawdę działa znacznie szybciej niż wersje svd w MATLAB. Możesz go pobrać tutaj:

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

Aby to przetestować, po prostu utwórz obraz w tym samym folderze (podobnie jak duża matryca, możesz sam stworzyć matrycę)

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

Szybka SVD

Kiedy uruchamiam go na pulpicie dla obrazu o rozmiarze 635 * 483, otrzymuję

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

Jak widać, przy niskich wartościach kjest on ponad 10 razy szybszy niż przy użyciu Matlab SVD. Nawiasem mówiąc, może być potrzebna następująca prosta funkcja dla funkcji testowej:

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

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

Nie dodałem metody PCA, ponieważ można ją łatwo wdrożyć za pomocą SVD. Możesz sprawdzić ten link, aby zobaczyć ich związek.

petrichor
źródło
12

możesz spróbować użyć kilku opcji.

1- Rozkład kar z matrycą . Stosujesz pewne ograniczenia karne wobec u i v, aby uzyskać odrobinę rzadkości. Szybki algorytm zastosowany w danych genomicznych

Zobacz Whitten Tibshirani. Mają także R-pkg. „Karany rozkład macierzy z aplikacjami do rzadkich głównych komponentów i kanoniczną analizą korelacji”.

2- Randomizowane SVD . Ponieważ SVD jest algorytmem nadrzędnym, pożądane może być bardzo szybkie przybliżenie, szczególnie w przypadku analizy eksploracyjnej. Za pomocą randomizowanego SVD możesz wykonywać PCA na ogromnych zestawach danych.

Zobacz Martinsson, Rokhlin i Tygert „Losowy algorytm rozkładu macierzy”. Tygert ma kod do bardzo szybkiej implementacji PCA.

Poniżej znajduje się prosta implementacja randomizowanego SVD w 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
źródło
+1 za karany rozkład macierzy. Ten pakiet jest niesamowity. Powinienem chyba wspomnieć, że jest napisane „Witten”, na wypadek, gdyby ludzie mieli problemy ze znalezieniem cytatu. Wreszcie, OP powiedział, że nie chcą niczego napisanego w R, ale zasadniczo każdy duży pakiet SVD będzie miał backend C, C ++ lub Fortran dla szybkości.
David J. Harris,
4

Wygląda na to, że chcesz użyć algorytmu Lanczos . W przeciwnym razie możesz skonsultować się z Golubem i Van Loanem. Kiedyś kodowałem algorytm SVD (w SML wszystkich języków) z ich tekstu i działał on całkiem dobrze.

shabbychef
źródło
3

Sugeruję wypróbowanie jądra PCA, które ma złożoność czasowo-przestrzenną zależną od liczby przykładów (N), a nie od liczby funkcji (P), które moim zdaniem byłyby bardziej odpowiednie w twoim ustawieniu (P >> N)). Jądro PCA działa w zasadzie z macierzą jądra NxN (macierz podobieństw między punktami danych), a nie z macierzą kowariancji PxP, z którą trudno jest sobie poradzić w przypadku dużych P. Kolejną dobrą rzeczą w jądrze PCA jest to, że może nauczyć się projekcji nieliniowych również jeśli używasz go z odpowiednim jądrem. Zobacz ten artykuł na temat PCA jądra .

ebony1
źródło
2

Wydaje mi się, że można wykonać PCA, obliczając rozkład własny X ^ TX zamiast XX ^ T, a następnie przekształcić, aby uzyskać komputery PC. Jednak nie pamiętam szczegółów z ręki, ale są one w (doskonałej) książce Jolliffe i sprawdzę je, kiedy będę w pracy. Transliterowałbym procedury algebry liniowej z np. Metod numerycznych w C, zamiast używać jakiegokolwiek innego algorytmu.

Dikran Torbacz
źródło
5
Dobry żal ... konstruowanie macierzy kowariancji nigdy nie jest najlepszym sposobem na SVD. Pokazałem przykład, dlaczego jawne tworzenie macierzy kowariancji nie jest dobrym pomysłem na matematyce . SE: math.stackexchange.com/questions/3869/3871#3871 .
JM nie jest statystykiem
1

Istnieje również metoda ładowania początkowego autorstwa Fisher i in. , Zaprojektowana dla kilkuset próbek o dużych wymiarach.

Główną ideą tej metody jest sformułowanie: „ponowne próbkowanie jest transformacją niskiego wymiaru”. Tak więc, jeśli masz małą (kilkaset) liczbę wysokowymiarowych próbek, nie możesz uzyskać więcej głównych składników niż liczba twoich próbek. Dlatego sensowne jest rozważenie próbek jako podstawy oszczędnej, rzutowanie danych na podprzestrzeń liniową obejmującą te wektory i obliczenie PCA w tej mniejszej podprzestrzeni. Dostarczają również więcej szczegółów na temat postępowania w przypadku, gdy nie wszystkie próbki mogą być przechowywane w pamięci.

Kolya Ivankov
źródło
0

Zobacz artykuł Sama Roweisa, EM Algorytmy dla PCA i SPCA .

ars
źródło
Algorytm Wikipedii cytuje to i jest równoważny z tym w przypadku znalezienia jednego głównego elementu na raz.
dsimcha
OK, widzę teraz link. Jest to dość proste podejście i, jak wspomina Wikipedia, istnieją postępy w tej podstawowej idei. Po namyśle będziesz musiał poradzić sobie z pewnymi kompromisami (w tym przypadku konwergencją). Zastanawiam się, czy zadajesz tutaj właściwe pytanie. Czy naprawdę nie ma dobrych powiązań z bibliotekami linalg dla D?
ars