54
$\begingroup$

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.

$\endgroup$
  • 2
    $\begingroup$ There are loads of public SVD algorithms. See en.wikipedia.org/wiki/…. Can't you use or adapt one of them? Also, R is open-source, and under a GPL licence, so why not borrow its algorithm if it does the job? $\endgroup$ – Rob Hyndman Sep 18 '10 at 3:16
  • $\begingroup$ @Rob: I'd like to avoid practically writing a linear algebra library, and I also want to avoid the GPL's copyleft. Also, I've looked at bits and pieces of the R source code before and it's generally not very readable. $\endgroup$ – dsimcha Sep 18 '10 at 3:25
  • 4
    $\begingroup$ Am I missing something? You have >10K features but <1K samples? This means the last 9K components are arbitrary. Do you want all 1K of the first components? $\endgroup$ – shabbychef Sep 18 '10 at 3:59
  • 2
    $\begingroup$ In any event, you can't escape having to implement SVD, though thanks to much numerical linear algebra research, there are now a lot of methods to choose from, depending on how big/small, sparse/dense your matrix is, or if you want just the singular values, or the complete set of singular values and left/right singular vectors. The algorithms are not terribly hard to understand IMHO. $\endgroup$ – J. M. is not a statistician Sep 18 '10 at 5:39
  • $\begingroup$ Can you tell us why you want to do PCA ? $\endgroup$ – robin girard Feb 13 '11 at 19:54
27
$\begingroup$

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

Fast SVD

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.

$\endgroup$
11
$\begingroup$

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)
}
$\endgroup$
  • $\begingroup$ +1 for penalized matrix decomposition. That package is pretty amazing. I should probably mention that it's spelled "Witten," though, in case people have trouble finding the cite. Lastly, the OP said they didn't want anything written in R, but essentially any large SVD package out there will have a C, C++, or Fortran backend for speed. $\endgroup$ – David J. Harris Jun 15 '11 at 13:41
4
$\begingroup$

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.

$\endgroup$
3
$\begingroup$

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.

$\endgroup$
2
$\begingroup$

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.

$\endgroup$
1
$\begingroup$

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.

$\endgroup$

We're looking for long answers that provide some explanation and context. Don't just give a one-line answer; explain why your answer is right, ideally with citations. Answers that don't include explanations may be removed.

0
$\begingroup$

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

$\endgroup$
  • $\begingroup$ The Wikipedia algorithm cites this and is equivalent to this for the case of finding one principal component at a time. $\endgroup$ – dsimcha Sep 18 '10 at 3:26
  • $\begingroup$ OK, I see the link now. This is a fairly simple approach, and like Wikipedia mentions, there are advances upon this basic idea. On reflection though, you're going to have to deal with some kind of trade-offs (convergence in this case). I wonder if you're asking the right question here. Are there really no good bindings to linalg libraries for D? $\endgroup$ – ars Sep 18 '10 at 20:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.