¿El mejor algoritmo PCA para un gran número de características (>10K)?

Question

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 la estadística y la programación.

Necesito escribir un código para hacer PCA (Análisis de Componentes Principales). He navegado por los algoritmos conocidos y he implementado este que, por lo que sé, es equivalente al algoritmo NIPALS. Funciona bien para encontrar los primeros 2-3 componentes principales, pero luego parece volverse 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 de un gran número de características (del orden de 10.000 a 20.000) y de tamaños de muestra del orden de unos pocos cientos.

2. Debe ser razonablemente implementable sin una biblioteca de álgebra lineal/matriz decente, ya que el lenguaje de destino es D, que aún no tiene una, e incluso si la tuviera, preferiría no añadirla como 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ápidamente, pero utiliza la descomposición del valor singular, que no es algo que quiera codificar yo mismo.

Answer 1

He implementado el SVD aleatorio como se indica en "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. Recuperado el 1 de abril de 2011, de http://arxiv.org/abs/1007.5510 .". Si quieres obtener el SVD truncado, realmente funciona mucho más rápido que las variaciones del svd en MATLAB. Puedes conseguirlo aquí:

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, basta con crear una imagen en la misma carpeta (al igual que 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 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 se puede ver, para valores bajos de k es más de 10 veces más rápido que utilizando el SVD de Matlab. 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 he añadido el método PCA, ya que es fácil de implementar utilizando SVD. Usted puede consulte este enlace para ver su relación.

Answer 2

Podrías intentar usar un par de opciones.

1- Descomposición matricial penalizada . Se aplican algunas restricciones de penalización en las u's y v's para obtener algo de sparsity. Algoritmo rápido que se ha utilizado en datos genómicos

Véase Whitten Tibshirani. También tienen un paquete R. "Una descomposición matricial penalizada, con aplicaciones a componentes principales dispersos y análisis de correlación canónica".

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

Véase Martinsson, Rokhlin y Tygert "A randomized algorithm for the decomposition of matrices". Tygert tiene código para una implementación muy rápida de PCA.

A continuación se muestra una implementación sencilla de la SVD aleatoria 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)
}

Answer 3

Parece que tal vez quieras usar el Algoritmo de Lanczos . En su defecto, podría consultar Golub & Van Loan. Una vez codifiqué un algoritmo SVD (en SML, de todos los lenguajes) a partir de su texto, y funcionó razonablemente bien.

Answer 4

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

Answer 5

Creo recordar que es posible realizar el ACP calculando la descomposición propia de X^TX en lugar de XX^T y luego transformar para obtener los PC. Sin embargo, no puedo recordar los detalles, pero está en el (excelente) libro de Jolliffe y lo buscaré la próxima vez que esté en el trabajo. Yo transliteraría las rutinas de álgebra lineal de, por ejemplo, Numerical Methods in C, en lugar de utilizar cualquier otro algoritmo.