# How to perform PCA for data of very high dimensionality?

To perform principal component analysis (PCA), you have to subtract the means of each column from the data, compute the correlation coefficient matrix and then find the eigenvectors and eigenvalues. Well, rather, this is what I did to implement it in Python, except it only works with small matrices because the method to find the correlation coefficient matrix (corrcoef) doesn't let me use an array with high dimensionality. Since I have to use it for images, my current implementation doesn't really help me.

I've read that it's possible to just take your data matrix $D$ and compute $DD^\top/n$ instead of $D^\top D/n$, but that doesn't work for me. Well, I'm not exactly sure I understand what it means, besides the fact that it's supposed to be a $n \times n$ matrix instead of $p\times p$ (in my case $p\gg n$). I read up about those in the eigenfaces tutorials but none of them seemed to explain it in such a way I could really get it.

In short, is there a simple algorithmic description of this method so that I can follow it?

• What you read is correct. The matrix $DD^\top$ is called the Gram matrix. Its eigenvectors are (scaled) principal components. Its eigenvalues are exactly identical, up to the factor $1/n$, to the eigenvalues of the covariance matrix $D^\top D/n$. Commented Feb 7, 2015 at 22:18

The easiest way to do standard PCA is to center the columns of your data matrix (assuming the columns correspond to different variables) by subtracting the column means, and then perform an SVD. The left singular vectors, multiplied by the corresponding singular value, correspond to the (estimated) principal components. The right singular vectors correspond to the (estimated) principal component directions — these are the same as the eigenvectors given by PCA. The singular values correspond to the standard deviations of the principal components (multiplied by a factor of root n, where n is the number of rows in your data matrix) — the same as the square root of the eigenvalues given by PCA.

If you want to do PCA on the correlation matrix, you will need to standardize the columns of your data matrix before applying the SVD. This amounts to subtracting the means (centering) and then dividing by the standard deviations (scaling).

This will be the most efficient approach if you want the full PCA. You can verify with some algebra that this gives you the same answer as doing the spectral decomposition of the sample covariance matrix.

There are also efficient methods for computing a partial SVD, when you only need a few of the PCs. Some of these are variants of the power iteration. The Lanczos algorithm is one example that is also related to partial least squares. If your matrix is huge, you may be better off with an approximate method. There are also statistical reasons for regularizing PCA when this is the case.

• Correct me if I am wrong, but I think Lanczos algorithm performs eigendecomposition and not SVD. Commented Feb 7, 2015 at 22:11
• An interested reader can look here for the further details on performing PCA via SVD: Relationship between SVD and PCA. How to use SVD to perform PCA? Commented Feb 7, 2015 at 22:14
• Is there any reason to not divide with standard deviations (scaling) prior to using PCA? Commented Jan 3 at 13:40

What you're doing right now is close, but you need to make sure you multiply the eigenvectors of (data . data.T) / lines on the left by data.T, in order to get the eigenvectors of (data.T . data) / lines. This is sometimes called the "transpose trick".

Here are some more details. Suppose you have a matrix $A$ that you want to perform PCA on; for simplicity, suppose that the columns of $A$ have already been normalized to have zero mean, so that we just need to compute the eigenvectors of the covariance matrix $A^T A$.

Now if $A$ is an $m \times n$ matrix, with $n >> m$, then $A^T A$ is a very large $n \times n$ matrix. So instead of computing the eigenvectors of $A^T A$, we might like to compute the eigenvectors of the much smaller $m \times m$ matrix $A A^T$ -- assuming we can figure out a relationship between the two. So how are the eigenvectors of $A^T A$ related to the eigenvectors of $A A^T$?

Let $v$ be an eigenvector of $A A^T$ with eigenvalue $\lambda$. Then

• $AA^T v = \lambda v$
• $A^T(A A^T v) = A^T(\lambda v)$
• $(A^T A)(A^T v) = \lambda (A^T v)$

In other words, if $v$ is an eigenvector of $A A^T$, then $A^T v$ is an eigenvector of $A^T A$, with the same eigenvalue. So when performing a PCA on $A$, instead of directly finding the eigenvectors of $A^T A$ (which may be very expensive), it's easier to find the eigenvectors $v$ of $AA^T$ and then multiply these on the left by $A^T$ to get the eigenvectors $A^T v$ of $A^T A$.

• This sounds like the "kernel trick" applied to PCA. en.wikipedia.org/wiki/Kernel_PCA It's a very good way of handling certain large matrices. Commented Feb 14, 2011 at 23:11
• +1. Perhaps one should add that $AA^\top$ is called the Gram matrix. Commented Feb 7, 2015 at 22:16
• To add on the answer of @raegtin, it should be noted that if $v$ is a normalized eigenvector of $A A^\top$, $A^\top v$ is an unormalized eigenvector of $A^\top A$. Commented Mar 26 at 16:55

It sounds like what you want is the NIPALS algorithm for performing PCA. It's a very popular algorithm among statisticians. It has many advantages:

• Computationally less expensive than SVD or eigenvalue decomposition methods if only the first few components are required.
• Has more modest storage requirements in general because the covariance matrix is never formed. This is a very important property for very large datasets.
• Can handle missing data in the dataset (though that's not an issue in your problem, since you're dealing with images).

Algorithm
Here's a simple and excellent description of the algorithm (in section 1.2)
http://stats4.eng.mcmaster.ca/w/mediafiles/mediawiki/f/f7/Section-Extra-Class-1.pdf

Remember to mean-center-scale first before doing PCA as it is scale-sensitive.

To add on Gilead's answer, they are computationally less expensive algorithms for truncated PCAs. NIPALS is indeed very popular, but I have had a lot of success with approximate methods that perform a succession of fits on partial data (what is often called PCA by random projection). This was discussed in a metaoptimize thread.

As you mention Python, let me point out that the algorithm is implemented in the scikit-learn: the PCA class. In particular, it is used in an example demonstrating eigenfaces.

If it helps someone, just adding on top of @raegtin's post, the following R code demonstrates how high-dimensional PCA can be done faster and more efficiently.

For the following matrix $$A_{m\times n}$$ with $$n>>m$$, the matrix $$A^TA$$ is a $$n\times n$$ matrix and not a full-rank matrix, having rank $$=m$$, so that only the first $$m$$ eigenvalues are non-zero.

Whereas $$AA^T$$ is a $$m \times m$$ matrix, which is a full-rank matrix, again with rank $$m$$ and $$A^Tv$$ is an eigenvector of $$A^TA$$ when $$v$$ is an eigenvector of $$AA^T$$. But we have to keep in mind that $$A^Tv$$ computed that way will not be of unit length, so we need to normalize the eigenvectors computed this way.

m <- 10
n <- 1000
A <- matrix(rnorm(m*n), m, n)

AtA <- t(A) %*% A
dim(AtA)
# [1] 1000 1000
qr(AtA)$rank # [1] 10 AAt <- A %*% t(A) dim(AAt) # [1] 10 10 qr(AAt)$rank
# [1] 10

# eigenvalues and eigenvectors
res <- eigen(AtA)
val1 <- res$$values vec1 <- res$$vectors

res <- eigen(AAt)
val2 <- res$$values vec2 <- res$$vectors
vec2 <- t(A) %*% vec2

# eigenvalues
length(val1)
# [1] 1000
length(val2)
# [1] 10
val1 <- round(val1, 3)
val2 <- round(val2, 3)
val1
# [1] 1154.743 1098.774 1075.320 1049.960 1035.392  996.713  943.589  924.477 886.510  814.518
# 0.000    0.000    0.000    0.000    0.000   0.000 ...
val2
# [1] 1154.743 1098.774 1075.320 1049.960 1035.392  996.713  943.589  924.477  886.510  814.518

all.equal(val1[1:10], val2)
# [1] TRUE

# eigenvectors
sqrt(sum(vec1[,1]^2))
# [1] 1
sqrt(sum(vec2[,1]^2))
# [1] 33.9815

# normalize eigenvectors
vec2 <- apply(vec2, 2, function(x) x / sqrt(sum(x^2)))
sqrt(sum(vec2[,1]^2))
# [1] 1

# compare run times
library(microbenchmark)
library(ggplot2)

mbm <- microbenchmark("PCA_AtA" = {
res <- eigen(AtA)
val1 <- res$$values vec1 <- res$$vectors
},
"PCA_AAt" = {
res <- eigen(AAt)
val2 <- res$$values vec2 <- res$$vectors
vec2 <- t(A) %*% vec2
})

autoplot(mbm)