¿El mejor algoritmo de PCA para una gran cantidad de características (> 10K)?

Anteriormente pregunté esto en StackOverflow, pero parece que podría ser más apropiado aquí, dado que no obtuvo ninguna respuesta en SO. Es una especie de intersección entre estadísticas y programación.

Necesito escribir un código para hacer PCA (Análisis de componentes principales). He examinado los algoritmos conocidos e implementado este , que por lo que puedo decir es equivalente al algoritmo NIPALS. Funciona bien para encontrar los primeros 2-3 componentes principales, pero luego parece ser muy lento para converger (del orden de cientos a miles de iteraciones). Aquí están los detalles de lo que necesito:

1. El algoritmo debe ser eficiente cuando se trata con un gran número de características (orden de 10,000 a 20,000) y tamaños de muestra del orden de unos pocos cientos.

2. Debe ser razonablemente implementable sin una biblioteca de álgebra / matriz lineal decente, ya que el lenguaje de destino es D, que aún no tiene uno, e incluso si lo tuviera, preferiría no agregarlo como una dependencia al proyecto en cuestión .

Como nota al margen, en el mismo conjunto de datos, R parece encontrar todos los componentes principales muy rápido, pero utiliza una descomposición de valores singulares, que no es algo que quiera codificar yo mismo.

Implementé la SVD aleatoria como se indica en "Halko, N., Martinsson, PG, Shkolnisky, Y., y Tygert, M. (2010). Un algoritmo para el análisis de componentes principales de grandes conjuntos de datos. Arxiv preprint arXiv: 1007.5510, 0526. Recuperado el 1 de abril de 2011, de http://arxiv.org/abs/1007.5510 ". Si desea obtener SVD truncado, realmente funciona mucho más rápido que las variaciones de svd en MATLAB. Puedes obtenerlo aqui:

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

Para probarlo, simplemente cree una imagen en la misma carpeta (como una matriz grande, puede crear la matriz usted mismo)

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

Cuando lo ejecuto en mi escritorio para obtener una imagen de tamaño 635 * 483, obtengo

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

Como puede ver, para valores bajos de k, es más de 10 veces más rápido que usar Matlab SVD. Por cierto, es posible que necesite la siguiente función simple para la función de prueba:

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

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

No agregué el método PCA ya que es fácil de implementar usando SVD. Puede consultar este enlace para ver su relación.

podrías intentar usar un par de opciones.

1- Descomposición matricial penalizada . Aplica algunas restricciones de penalización en las u y las v para obtener cierta dispersión. Algoritmo rápido que se ha utilizado en datos genómicos.

Ver Whitten Tibshirani. También tienen un R-pkg. "Una descomposición de matriz penalizada, con aplicaciones para dispersar componentes principales y análisis de correlación canónica".

2- SVD aleatorizado . Dado que SVD es un algoritmo maestro, puede ser conveniente encontrar una aproximación muy rápida, especialmente para el análisis exploratorio. Usando SVD aleatorio, puede hacer PCA en grandes conjuntos de datos.

Ver Martinsson, Rokhlin y Tygert "Un algoritmo aleatorio para la descomposición de matrices". Tygert tiene código para una implementación muy rápida de PCA.

A continuación se muestra una implementación simple de SVD aleatorizada en 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)
}

Parece que tal vez quieras usar el algoritmo de Lanczos . De lo contrario, es posible que desee consultar Golub & Van Loan. Una vez codifiqué un algoritmo SVD (en SML, de todos los idiomas) de su texto, y funcionó razonablemente bien.

Sugeriría probar el kernel PCA que tiene una complejidad de tiempo / espacio que depende de la cantidad de ejemplos (N) en lugar de la cantidad de características (P), lo que creo que sería más adecuado en su entorno (P >> N)). Kernel PCA básicamente funciona con la matriz de kernel NxN (matriz de similitudes entre los puntos de datos), en lugar de la matriz de covarianza PxP, que puede ser difícil de manejar para P. grande. Otra cosa buena sobre la PCA de kernel es que puede aprender proyecciones no lineales también si lo usa con un núcleo adecuado. Vea este documento sobre el kernel PCA .

Me parece recordar que es posible realizar PCA calculando la descomposición propia de X ^ TX en lugar de XX ^ T y luego transformar para obtener las PC. Sin embargo, no puedo recordar los detalles de antemano, pero está en el libro (excelente) de Jolliffe y lo buscaré la próxima vez que esté en el trabajo. Transliteraría las rutinas de álgebra lineal de, por ejemplo, métodos numéricos en C, en lugar de usar cualquier otro algoritmo.

También existe el método bootstrap de Fisher et al , diseñado para varios cientos de muestras de alta dimensión.

La idea principal del método se formula como "el remuestreo es una transformación de baja dimensión". Entonces, si tiene un número pequeño (varios cientos) de muestras de alta dimensión, entonces no puede obtener más componentes principales que el número de sus muestras. Por lo tanto, tiene sentido considerar las muestras como una base parsimoniosa, proyectar los datos en el subespacio lineal que abarcan estos vectores y calcular PCA dentro de este subespacio más pequeño. También proporcionan más detalles sobre cómo tratar el caso cuando no todas las muestras se pueden almacenar en la memoria.

Ver el artículo de Sam Roweis, Algoritmos EM para PCA y SPCA .