Best PCA algorithm for huge number of features (>10K)?

Question

I previously asked this on StackOverflow, but it seems like it might be more appropriate here, given that it didn't get any answers on SO. It's kind of at the intersection between statistics and programming.

I need to write some code to do PCA (Principal Component Analysis). I've browsed through the well-known algorithms and implemented this one, which as far as I can tell is equivalent to the NIPALS algorithm. It works well for finding the first 2-3 principal components, but then seems to become very slow to converge (on the order of hundreds to thousands of iterations). Here are the details of what I need:

1. The algorithm must be efficient when dealing with huge numbers of features (order 10,000 to 20,000) and sample sizes on the order of a few hundred.

2. It must be reasonably implementable without a decent linear algebra/matrix library, as the target language is D, which doesn't have one yet, and even if it did, I would prefer not to add it as a dependency to the project in question.

As a side note, on the same dataset R seems to find all principal components very fast, but it uses singular value decomposition, which is not something I want to code myself.

Answer 1

I've implemented the Randomized SVD as given in "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.". If you want to get truncated SVD, it really works much much faster than the svd variations in MATLAB. You can get it here:

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

To test it, just create an image in the same folder (just as a big matrix,you can create the matrix yourself)

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

When I run it on my desktop for an image of size 635*483, I get

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

As you can see, for low values of k, it is more than 10 times faster than using Matlab SVD. By the way, you may need the following simple function for the test function:

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

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

I didn't add the PCA method since it is straightforward to implement using SVD. You may check this link to see their relationship.

Answer 2

you could trying using a couple of options.

1- Penalized Matrix Decomposition. You apply some penalty constraints on the u's and v's to get some sparsity. Quick algorithm that has been used on genomics data

See Whitten Tibshirani. They also have an R-pkg. " A penalized matrix decomposition, with applications to sparse principal components and canonical correlation analysis."

2- Randomized SVD. Since SVD is a master algorithm, find a very quick approximation might be desirable, especially for exploratory analysis. Using randomized SVD, you can do PCA on huge datasets.

See Martinsson, Rokhlin, and Tygert "A randomized algorithm for the decomposition of matrices". Tygert has code for a very fast implementation of PCA.

Below is a simple implementation of randomized SVD in 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

It sounds like maybe you want to use the Lanczos Algorithm. Failing that, you might want to consult Golub & Van Loan. I once coded a SVD algorithm (in SML, of all languages) from their text, and it worked reasonably well.

Answer 4

I'd suggest trying kernel PCA which has a time/space complexity dependent on the number of examples (N) rather than number of features (P), which I think would be more suitable in your setting (P>>N)). Kernel PCA basically works with NxN kernel matrix (matrix of similarities between the data points), rather than the PxP covariance matrix which can be hard to deal with for large P. Another good thing about kernel PCA is that it can learn non-linear projections as well if you use it with a suitable kernel. See this paper on kernel PCA.

Answer 5

I seem to recall that it is possible to perform PCA by computing the eigen-decomposition of X^TX rather than XX^T and then transform to get the PCs. However I can't remember the details off-hand, but it is in Jolliffe's (excellent) book and I'll look it up when I am next at work. I'd transliterate the linear algebra routines from e.g. Numerical Methods in C, rather than use any other algorithm.

Answer 6

See Sam Roweis' paper, EM Algorithms for PCA and SPCA.

Answer 7

There is also the bootstrap method by Fisher et al, designed for several hundred samples of high dimension.

The main idea of the method is formulated as "resampling is a low-dimension transformation". So, if you have a small (several hundred) number of high-dimensional samples, then you can't get more principal components than the number of your samples. It thus makes sense to consider the samples as a parsimonious basis, project the data on the linear subspace spanned by these vectors, and calculate PCA within this smaller subspace. They also provide more details how to deal with the case when not all samples may be stored in the memory.