Principal component analysis

Joshua R. Davis, Carleton College, Math 134, Fall 2024

Introduction

This first chunk of code is used to make point clouds of data. (Activate the cell of code, but you are not supposed to understand how it works.)

Clear[cloud];
cloud[mu_, sigmaSq_, n_] :=
Table[
RandomVariate[
MultinormalDistribution[
mu, 
sigmaSq IdentityMatrix[Length[mu]]]], 
{n}];

Here is a data set consisting of 200 points in 3D. Upon first glance, it seems that the data are unimodal, meaning that they form one region of high density. But if you rotate the plot, you can see that the data are bimodal, meaning that they form two distinct high-density regions.

first = cloud[{0.8, 0, 1}, 0.01, 100];
second = cloud[{0.2, 1, 0}, 0.01, 100];
data = Join[first, second];
ListPointPlot3D[data, BoxRatios->{1, 1, 1}]

Statisticians care very much about this sort of thing. For starters, the mean of a bimodal data set tends to lie in the empty space between the two modes. In that case, the mean is not really a good summary of the data set, is it?

Here is a data set consisting of 200 points in 4D. To picture it, I project it down to three dimensions (the first three). The data look unimodal, even if you rotate the 3D plot.

first = cloud[{0, 0, 0, -1}, 0.01, 100];
second = cloud[{0, 0, 0, 1}, 0.01, 100];
data = Join[first, second];
dataProjected = data[[All, 1;;3]];
ListPointPlot3D[dataProjected, BoxRatios->{1, 1, 1}]

But look what happens if we instead project down the last three dimensions. Two modes are revealed.

dataProjected = data[[All, 2;;4]];
ListPointPlot3D[dataProjected, BoxRatios->{1, 1, 1}]

Many scientific data sets take the form of point clouds in high dimensions. For example, an ecologist friend of mine once asked me to analyze a 28-dimensional data set. The 28 dimensions arose, because she was measuring 28 distinct markers of paleoclimate. We can’t visualize 28 dimensions, but we still want to use our visual intuition to guide our study of such data. We have to project down to two or three dimensions, to make any kind of plot.

So the question is: Given a high-dimensional data set, what is the best “angle” from which to view it, so that its structure is manifested most clearly?

Sample mean and covariance

Here is some code for computing the sample mean and sample covariance matrix of a data set. (As always, you should execute the code cell, so that you can use it, but you aren’t expected to understand the code itself.)

Clear[mean, covariance];
mean[data_] :=
Total[data] / Length[data];
covariance[data_] :=
Module[{m, xHats, b}, 
m = mean[data];
xHats = Map[Function[x, x - m], data];
b = Transpose[xHats];
b.Transpose[b] / (Length[data] - 1)];

Assuming that you still have the 4D data set from above in Mathematica’s memory (if not, then re-generate it), here is its 4D sample mean and 4x4 sample covariance.

MatrixForm[mean[data]]
MatrixForm[covariance[data]]

Of course, Mathematica has built-in functions to compute this stuff, because they’re so prevalent. I wrote my own versions above, just because I’m trying to de-mystify the whole process for you.

MatrixForm[Mean[data]]
MatrixForm[Covariance[data]]

Principal component analysis (PCA)

Principal component analysis (PCA) might sound scary, but essentially it means: Diagonalize the sample covariance.

For the sample covariance just computed, here are the eigenvalues and the eigenvectors as columns in a matrix P. Here’s what they mean. The eigenvalues measure how spread out the data set is, and the eigenvectors describe the directions of that spreading.

{eigvals, eigvecs} = Eigensystem[Covariance[data]];
eigvals
p = Transpose[eigvecs];
MatrixForm[p]

What is the following code checking? And can you guess what Chop means?

MatrixForm[Transpose[p].p]
MatrixForm[Chop[Transpose[p].p]]

To get back to our original motivating question --- how to plot high-dimensional data --- first we shift the data so that they have mean zero.

shifted = Map[Function[x, x - Mean[data]], data];
Mean[shifted]
Chop[Mean[shifted]]

Then we rotate (and possibly flip) those shifted data using the orthogonal matrix P^T. Then we project onto the first two dimensions and plot. The two modes are apparent. If you examine the tick marks, then you see that the data are more spread in the horizontal direction than in the vertical direction.

rotated = Map[Function[x, Transpose[p].x], shifted];
ListPlot[rotated[[All, 1;;2]], AspectRatio-> 1]

Here’s what’s happened. The sample covariance is telling us what the biggest directions of spreading are. The plot above is showing just the two biggest directions of spreading. It does this by rotating (after shifting and possibly flipping) the data, before projecting them down to two dimensions. In essence, it is viewing the data from the best possible viewpoint.

If you were really hoping to end up with a 3D plot, then we can do that instead. It’s essentially just tagging one more dimension onto the 2D plot above.

ListPointPlot3D[rotated[[All, 1;;3]], BoxRatios->{1, 1, 1}]

SVD?

What does any of this have to do with the SVD? Well, remember that PCA is all about diagonalizing the sample covariance S = (B B^T) / (n - 1). If we let A = B^T / sqrt(n - 1), then PCA is all about diagonalizing A^T A. But diagonalizing A^T A is essentially the same thing as computing the SVD of A, and the SVD approach is more numerically robust. So, in actual statistics software, PCA is computed through the SVD (as far as I understand).

To drive home this point, here are the same PCA computations, done through the SVD, without even computing the sample covariance.

shifted = Map[Function[x, x - Mean[data]], data];
a = Transpose[shifted] / Sqrt[Length[data] - 1];
{u, s, v} = SingularValueDecomposition[a];
rotated =  Map[Function[x, Transpose[u].x], shifted];
ListPlot[rotated[[All, 1;;2]], AspectRatio-> 1]

Clustering

What else can we do with data? Lots of stuff. For example, because these data are bimodal, maybe we’d like to run a clustering algorithm, to split the data into two groups. (How this works, and the many choices that are being made behind the scenes, is a big topic. Let’s not go into it.)

clusters = GatherBy[data, ClusterClassify[data]];
Length[clusters]

I hope that you got two clusters above. (Rarely, the results are bad. In that case, try starting over with a new data set.) Here is a naive plot of the two clusters --- one in red and one in blue --- showing the first three dimensions. From this viewpoint, the choice of what’s red and what’s blue doesn’t make much sense.

firstClus = clusters[[1]];
secondClus = clusters[[2]];
Show[
ListPointPlot3D[firstClus[[All, 1;;3]], BoxRatios->{1, 1, 1}, PlotStyle->Red], 
ListPointPlot3D[secondClus[[All, 1;;3]], BoxRatios->{1, 1, 1}, PlotStyle->Blue]]

But if we do the shifting, rotation, and maybe flipping of PCA, then the coloration makes lots of sense.

shifted = Map[Function[x, x - Mean[data]], data];
a = Transpose[shifted] / Sqrt[Length[data] - 1];
{u, s, v} = SingularValueDecomposition[a];
firstTransf = Map[Function[x, Transpose[u].(x - Mean[data])], firstClus];
secondTransf = Map[Function[x, Transpose[u].(x - Mean[data])], secondClus];
plotRange = {
MinMax[Join[firstTransf, secondTransf][[All, 1]]], 
MinMax[Join[firstTransf, secondTransf][[All, 2]]]};
Show[
ListPlot[firstTransf[[All, 1;;2]], AspectRatio-> 1, PlotStyle->Red, PlotRange->plotRange], 
ListPlot[secondTransf[[All, 1;;2]], AspectRatio-> 1, PlotStyle->Blue]]

Lossy data compression

Data compression is the idea of taking a chunk of data and finding a smaller representation for it, so that it uses less space. Compression is important in many computer applications. For example, when you stream video over the Internet from Netflix, the video is stored in a compressed format on Netflix’s servers. Netflix sends it to you compressed. Your video player decompresses it just before showing it to you.

No one compression scheme is appropriate for all applications. So there are many compression schemes, which fall into two big families. Lossless compression is perfect, in that the decompressed version exactly matches the pre-compression version. An example is ZIP, which is used to compress software applications, which must be perfect in order to function correctly. Lossy compression is approximate, in that the decompressed version is close to the pre-compression version, where “close to” depends on the application. For example, the JPEG compression system for images is lossy. For another example, streaming video is lossily compressed.

PCA gives us a simple way to do lossy data compression. It is appropriate for data where some of the dimensions can be ignored because the sample covariance is small in those dimensions.

Let’s do an example using randomly generated data. To get started, we pick the eigenvalues of the covariance. There are a few big ones and many small ones. How many there are total is the dimension d. And pick any large sample size n you like.

sigmaSqs = {0.99, 0.98, 0.97, 0.06, 0.05, 0.04, 0.03, 0.02, 0.01};
d = Length[sigmaSqs];
n = 1000;

This next chunk of code makes a covariance matrix, with the eigenvalues you prescribed above, but with a random orientation.

symm = Table[RandomVariate[NormalDistribution[]], {d} ,{d}];
symm = symm + Transpose[symm];
orthog = Eigenvectors[symm];
orthog = Map[Function[v, v / Norm[v]], orthog];
sigmaTrue = Transpose[orthog].DiagonalMatrix[sigmaSqs].orthog;

This next chunk of code picks a random mean for the data set. Then it uses that mean and covariance matrix to make a data set of n d-dimensional vectors.

mTrue = Table[RandomVariate[NormalDistribution[]], {d}];

Then, from the mean and covariance just chosen, we generate a data set.

data = Table[RandomVariate[MultinormalDistribution[mTrue, sigmaTrue]], {n}];

Here is the sample mean of the data set. Notice that it’s close to, but not exactly equal to, the mean that was used to construct the data.

mData = Mean[data];
MatrixForm[mData]
MatrixForm[mTrue]

Here is the sample covariance. Again, it is close to, but not exactly equal to, the covariance that was used to construct the data.

MatrixForm[Covariance[data]]
MatrixForm[sigmaTrue]

Now we do PCA. The matrix U tells us the principal directions, and the matrix S tells us the singular values.

shifted = Map[Function[x, x - mData], data];
a = Transpose[shifted] / Sqrt[Length[data] - 1];
{u, s, v} = SingularValueDecomposition[a];
MatrixForm[u]
MatrixForm[s]

The squares of the singular values are close to the covariance eigenvalues that we started with.

MatrixForm[s.Transpose[s]]
sigmaSqs

Now we compress the data set. For each datum x, we compute w = P^T (x - m) and store only the first few entries.

comps = Map[Function[shift, (Transpose[u].shift)[[1;;3]]], shifted];
Dimensions[comps]

Now let’s decompress the data set.

decomps = Map[Function[comp, u.Join[comp, Table[0, {6}]] + mData], comps];
Dimensions[decomps]

You can see that the first decompressed datum is close to, but not equal to, the first original datum. That’s lossy compression.

data[[1]]
decomps[[1]]