Best PCA algorithm for huge number of features (>10K)? 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:


*

*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.

*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.
 A: 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.
A: 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.
A: 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.
A: 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)
}

A: See Sam Roweis' paper, EM Algorithms for PCA and SPCA.
A: 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.
