Meilleur algorithme PCA pour un grand nombre de fonctionnalités (> 10K)?

J'ai déjà posé cette question à StackOverflow, mais il semble que cela conviendrait mieux ici, étant donné que cela n'a pas donné de réponse à SO. C'est un peu à la croisée des statistiques et de la programmation.

J'ai besoin d'écrire du code pour faire PCA (analyse en composantes principales). J'ai parcouru les algorithmes bien connus et mis en œuvre celui-ci , qui pour autant que je sache est équivalent à l'algorithme NIPALS. Cela fonctionne bien pour trouver les 2-3 premières composantes principales, mais semble ensuite devenir très lent à converger (de l'ordre de centaines à des milliers d'itérations). Voici les détails de ce dont j'ai besoin:

1. L'algorithme doit être efficace pour traiter un grand nombre de caractéristiques (ordre de 10 000 à 20 000) et des tailles d'échantillon de l'ordre de quelques centaines.

2. Il doit être raisonnablement implémentable sans une bibliothèque matricielle d'algèbre linéaire, car le langage cible est D, qui n'en a pas encore, et même s'il en possédait un, je préférerais ne pas l'ajouter comme dépendance du projet en question. .

En passant, sur le même jeu de données, R semble trouver toutes les composantes principales très rapidement, mais il utilise une décomposition en valeurs singulières, ce que je ne veux pas coder moi-même.

J'ai mis en œuvre la SVD randomisée comme indiqué dans "Halko, N., PG Martinsson, Y Shkolnisky, et M. Tygert (2010). Un algorithme pour l'analyse en composantes principales de grands ensembles de données. Arxiv preprint arXiv: 1007.5510, 0526. Extrait le 1 er avril 2011 de http://arxiv.org/abs/1007.5510 . ". Si vous voulez obtenir des fichiers SVD tronqués, cela fonctionne vraiment beaucoup plus rapidement que les variations de DVD dans MATLAB. Vous pouvez l'obtenir ici:

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

Pour le tester, créez simplement une image dans le même dossier (juste comme une grande matrice, vous pouvez créer la matrice vous-même)

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

Lorsque je l'exécute sur mon bureau pour une image de taille 635 * 483, je reçois

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

Comme vous pouvez le constater, pour les faibles valeurs de k, il est plus de 10 fois plus rapide que l’utilisation de Matlab SVD. En passant, vous aurez peut-être besoin de la fonction simple suivante pour la fonction de test:

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

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

Je n'ai pas ajouté la méthode PCA car il est facile de l'implémenter avec SVD. Vous pouvez vérifier ce lien pour voir leur relation.

vous pouvez essayer d'utiliser plusieurs options.

1- Décomposition de la matrice pénalisée . Vous appliquez des contraintes de pénalité sur les u et v pour obtenir un peu de parcimonie. Algorithme rapide utilisé pour les données génomiques

Voir Whitten Tibshirani. Ils ont aussi un R-pkg. "Une décomposition matricielle pénalisée, avec des applications pour les composantes principales éparses et l'analyse de corrélation canonique."

2- SVD randomisée . Comme la SVD est un algorithme maître, il peut être souhaitable de disposer d’une approximation très rapide, en particulier pour une analyse exploratoire. À l’aide de la SVD randomisée, vous pouvez effectuer une PCA sur d’énormes jeux de données.

Voir Martinsson, Rokhlin et Tygert "Un algorithme randomisé pour la décomposition de matrices". Tygert a un code pour une implémentation très rapide de la PCA.

Vous trouverez ci-dessous une implémentation simple de SVD randomisée dans 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)
}

Il semble que vous souhaitiez peut-être utiliser l’ algorithme Lanczos . Sinon, vous voudrez peut-être consulter Golub & Van Loan. Une fois, j'ai codé un algorithme SVD (en SML, dans toutes les langues) à partir de leur texte, et cela a fonctionné assez bien.

Je suggérerais d'essayer la PCA du noyau dont la complexité spatio-temporelle dépend du nombre d'exemples (N) plutôt que du nombre de fonctionnalités (P), ce qui, à mon avis, conviendrait mieux dans votre contexte (P >> N)). Le noyau PCA fonctionne essentiellement avec la matrice de noyau NxN (matrice de similitudes entre les points de données), plutôt que la matrice de covariance PxP, qui peut être difficile à gérer pour les grands P. Un autre avantage de la PCA noyau est qu'il peut apprendre des projections non linéaires. ainsi si vous l'utilisez avec un noyau approprié. Voir cet article sur la PCA du noyau .

Il me semble rappeler qu'il est possible d'effectuer une ACP en calculant la décomposition propre de X ^ TX plutôt que de XX ^ T, puis de la transformer pour obtenir les PC. Cependant, je ne me souviens pas des détails, mais c'est dans le livre (excellent) de Jolliffe et je le vérifierai quand je serai au travail. Je translittérerais les routines d'algèbre linéaire de Numericical Methods in C, par exemple, plutôt que d'utiliser un autre algorithme.

Il existe également la méthode bootstrap de Fisher et al , conçue pour plusieurs centaines d’échantillons de grande dimension.

L'idée principale de la méthode est formulée comme suit: "le rééchantillonnage est une transformation de faible dimension". Ainsi, si vous avez un petit nombre (plusieurs centaines) d'échantillons de grande dimension, vous ne pouvez pas obtenir plus de composants principaux que le nombre d'échantillons. Il est donc logique de considérer les échantillons comme une base parcimonieuse, de projeter les données sur le sous-espace linéaire couvert par ces vecteurs et de calculer la PCA dans ce sous-espace plus petit. Ils fournissent également plus de détails sur la manière de traiter le cas lorsque tous les échantillons ne peuvent pas être stockés dans la mémoire.

Voir l'article de Sam Roweis, EM Algorithms for PCA et SPCA .