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

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

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

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.

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)
}

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.

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 .

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.

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.

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