PCA / FA. A Map to my posts

The Posts

The purpose of this post is to provide guidance to a reader who has just discovered that I have a large pile of posts about principal components / factor analysis. This pile of posts might seem a very jungle, without any map.

Here, have a map.

As I finalize this post, it will be number 52 in PCA / FA. Here’s a list of the 52 posts, including the dates spanned by any group, and the number of posts in that group. (When the picture was taken, I didn’t know when this would be published. In fact, post 51 was scheduled but not yet published. Even more, post 51 did not even exist when the first picture was created.)

transition/attitude matrices is a post that is sometimes relevant when we discuss “new data” in PCA, but it is not in the PCA / FA category.

“tricky prepro” is short for “tricky preprocessing”, and discusses the combination of constant row sums and covariance or correlation matrix.

Posts which have no named author, such as “example 8″, have the tag, “Rip”.

Right now the tag cloud includes entries such as Harman and Davis, i.e. the authors of textbooks I used. This means you can see all the posts based on Harman (including those which include Davis) by clicking on the Harman tag. The number of tags shown, however, is limited to 45. In time, some or all of these author tags may disappear. What to do then?

This post shows you that the original Harman posts all occurred in January 2008, and that there were some posts about Harman & Davis in April 2008, and the “et al” on post 51 means that it talks about Harman, Jolliffe, and Basilevsky.. The dates would let you pick whatever you wanted out of the monthly archives.

The Authors

Let me summarize the sources quickly. There are three texts which specialize in PCA / FA in general: Harman, Jolliffe, Basilevsky. There is one text which specializes in PCA / FA in chemistry: Malinowski. There are three texts which specialize in data analysis in a particular field: chemistry, geology, social science (resp. Brereton, Davis, and Bartholomew et al.)

Harman’s primary focus is factor analysis, but I started with him because he had a simple example with data. He also has a lot of material on what they call rotations, but I did not look at it. If you need to understand this, I think he’s a good place to start. (What people mean in this field when they say “rotation” is anything but; in fact, they want to turn an orthogonal basis into a non-orthogonal basis, and that is not done by a rotation!)

Jolliffe is all about principal component analysis. As I said before, if you do this professionally, you probably need this book: he talks a lot about what you’re doing and what it means. I didn’t start here because his examples are not self-contained: he references the literature which contains the data.

Davis is about data analysis primarily for geology. His material on PCA / FA is limited to one chapter. Nevertheless, as you might guess from the number of posts, Davis is where I began to put it all together. I would recommend his book for data analysis in general, not just for PCA / FA.

Malinowski and Brereton are special-purpose texts, even within the category “PCA in chemistry” (to judge from the comment here.) From Malinowski I learned quite a bit more about the geometry of the SVD, the plausibility of using data which was not even centered, and an interesting approach to missing data.

Brereton has some nice data, but you have to own the book to access the data. Also, he did not look like the kind of book I could have learned this from. (Practice with, yes; learn the concepts, no. But then, I’m not trying to call a PCA command in a statistics package – not that it’s a bad thing to do, but it’s not my cup of tea.)

Bartholomew et al. is about data analysis primarily for the social sciences. Like Davis, their material on PCA is limited to one chapter, but with another on classical FA (which I have not discussed). The book, however, has publicly available data, and this makes it an excellent source for practice. In fact, it was here that I finally saw how easy it was to compute the scores. They do not use matrix notation, so it might help if you had learned that elsewhere before picking up this book.

Basilevsky is devoted to factor analytic methods: factor analysis, principal component analysis, and regression. He tends to have data printed in tables, so I have scanned things into the computer for myself in order to play with his examples. I think I’m just as glad I found him after I understood how to do this; he might have made me dizzy in the beginning.

Introduction to The Calculations

A summary of just about all the calculations I could do for PCA is here but you might also look at the post immediately preceding this one.

Now that I understand the basics of PCA, I suspect that I personally would use it in two ways. One, I might use it on data before I ran a regression. Two, I might do just what I have been doing, namely checking someone else’s calculations.

Using PCA Myself

Let’s take care of the simple case first. If I have some data, and if I choose to run PCA, I would be interested – at this stage of my understanding – in two things. One, how would a PCA redistribute the variance? Two, how many linearly independent variables do I have?

I would answer both of those questions by getting the SVD (singular value decomposition, ) of my data matrix X. The number of nonzero singular values w will tell me the rank of the matrix; the number of “small” singular values will let me estimate multicolinearity. Since I want to know about the redistributed variance, I need to use either centered or standardized data. (Otherwise, the eigenvalues are not equal to variances.) In either case, it is the orthogonal eigenvector matrix v (not the weighted matrix A) which will show me the redistributed variances.

When I do this, I’ll show it to you. (I have yet to do it, since understanding PCA.) Who knows? At that time, I may figure out that there are other things I would wish to get out of the PCA.

Anyway, it seems to me today that for my own purposes I would simply get the SVD of my data, and look at w and v. This is way, way simpler than what I’ve usually been doing while checking other peoples calculations.

(I do wonder if the matrix u can be used to isolate variables which are multi-collinear. We saw in Malinowski that we could construct a projection operator…. Let me think about that.)

The other case is somewhat more complicated. This takes us back to my beginnings.

If I am checking someone else’s analysis, the very first question is: did they give me the data?

Checking Someone Else, No Data

If they did not give me the data, they must have given me a “dispersion matrix”: some form of either or , in fact probably either a covariance matrix or a correlation matrix. Whatever it is, denote it by c.

I would then get an eigendecomposition of c,

,

where the eigenvalues of c are the nonzero elements of the diagonal matrix . To put that another way, the diagonal matrix has for its diagonal elements.

You will recall that minor hell can break loose. Although any matrix of the form or must be positive definite – when computed exactly – any rounded off representation of it need not be positive definite. We saw an example here. The specific problem we might encounter is that a very small positive eigenvalue – really – might come out negative when we compute using a printed – rounded-off – dispersion matrix instead of the computed matrix. (In principle, of course, even the computed one could fail to be positive definite; there’s nothing sacred about rounding to 6 places or 15 places instead of 3 or so.)

If that happens, I would set the troublesome eigenvalue to zero, and continue. Incidentally, this can occur if we started with constant-row-sum data and then computed either the covariance or correlation matrix. I discussed that here.

Moving on, let us suppose we have the orthogonal eigenvector matrix v and the diagonal matrix of nonnegative eigenvalues (i.e. we took care of any negative ones, by setting them to 0).

If the work we are checking stays with v and , then he is probably doing the kinds of things Jolliffe discussed.

Otherwise, I expect to see a weighted eigenvector matrix

.

If we were not given data, the only conclusions we can check are those drawn from v, A, and . Computing the “scores” requires the data, and even if the analysis shows the scores, we cannot check them.

But what if we were given the data? Well, this could get interesting. (This, after all, is why I’ve been doing this stuff for a year and a half!)

Checking Someone Else, With Data

The first significant consequence of being given the data is this: we should be able to relate their eigenvalues and eigenvectors to the data. (We’ll get to the scores….)

That may not be straightforward.

The first challenge is: are the variables in columns or in rows? Harman, we will recall, puts the variables in rows, using the transpose of our customary X matrix: hence we had his

and

.

I would compute the means and variances of each variable (i.e. of each column or row). If I were at all uncertain whether the variables were rows or columns, I would compute the means and variances of each row and column. In fact, even if I know that the variables are in columns, say, I would compute the row sums.

This tells me whether I was given raw data, centered data, standardized data, or small standard deviates; or something else. It remains to be seen if they preprocessed the data before computing the dispersion matrix.

Okay, if the given data is standardized (or small standard deviates), I expect that they used the given data for the analysis. But I might be wrong.

Now, what did they do? They either computed an SVD of the data, or an eigendecomposition of a dispersion matrix. I would compute whatever I think they did with my best guess for the preprocessing.

If my eigenvectors are right, then I have the correct dispersion matrix to within a scale factor; if my eigenvalues are also right, then the scale factor is right.

Bear in mind that any multiple of an eigenvector is also an eigenvector. I would have computed an orthogonal eigenvector matrix v. If the given eigenvectors are of unit length, we have a match so long as they agree with mine to within a sign. If the given eigenvectors do not agree with mine to within a sign, I would compute their lengths. If their lengths turn out to be , I would immediately compute

,

expecting that my columns of A would agree with theirs to within a sign. About the only other scaling I have seen is to set the largest component of each eigenvector to 1. (Apart from Basilevsky normalizing the rows of the A matrix by the inverse standard deviations! And that, by the way, means the columns of A are no longer eigenvectors.)

Let me remind you that raw data, centered data, and standardized data (each denoted X) lead to different eigenvectors of ; so if I match the eigenvectors, I have correctly matched the author’s raw, centered, or standardized data. Until I match the eigenvectors, there’s not much point in worrying about the eigenvalues.

If the eigenvalues are off by a scale factor, once the eigenvectors are right, then I need to use a different multiple of . That is,

, , and

have the same eigenvectors, but the eigenvalues differ by factors of N or N-1 or (N-1)/N.

The same considerations apply to the computation of the singular value decomposition. If my matrices u and v agree with the author (to within signs, assuming orthonormal eigenvectors, as usual), I have correctly matched raw, centered, or standardized. If the singular values w are off by a scale factor, then they scaled the data somehow.

Why all the fuss?

• I have seen people use 1/N in their theoretical discussion, but use 1/(N-1) in their computation.
• I have seen someone use 1/N to get “the correlation matrix” or “the covariance matrix”, but I would have used 1/(N-1).
• I have seen people say they were using “the covariance matrix”, but they actually used instead of .
• I have seen people call the covariance matrix when X wasn’t even centered.
• I have seen someone say the data was standardized, when in fact it was small standard deviates. (Same eigenvectors, different eigenvalues.)

Of course, either the author or I might have made a mistake getting the data into the computer; but I have learned to consider instead the possibility that an author and I are doing different computations on the same given data.

Finally, I have sometimes found that I agree with an author’s largest few eigenvectors (by which I mean, the eigenvectors associated with the largest eigenvalues), but not with the smaller ones. I suspect that this is a result of their using a sequential algorithm to find eigenvectors one after another.

I remind you that it is very easy to check an eigenvector-eigenvalue pair, my own or an author’s. Suppose we have done an eigendecomposition of a matrix c,

.

Suppose, for example, that is the third eigenvector, with associated eigenvalue ). Then the definitions of eigenvector and eigenvalue tell us that

.

The LHS is a matrix-vector product and the RHS is a scalar times a vector. Compute both sides and compare them.

Whew! So far, I have matched the author’s eigendecomposition. Presumably, I know whether they used the orthogonal eigenvector matrix v or the weighted eigenvector matrix A. The eigenvector matrix, whichever one they chose, is usually the “loadings”.

We’re actually home free, now.

If they present the “scores”, I expect to find that they are either X v or – if they come from A – the first few columns of ( here), unless I have a reason to believe they are using Davis’ scheme:

.

(And for more about that, see the post immediately preceding this one.)

PCA/FA Answers to some Basilevsky questions

Let us look at three of the questions I asked early in February, and answer two of them.

First, what do we know? What have we done?

We assume that we have data X, with the variables in columns, as usual. In fact, we assume that the data is at least centered, and possibly standardized.

We compute the covariance matrix

,

then its eigendecomposition

,

where is the diaginal matrix of eigenvalues. We define the -weighted matrix

.

Finally, we use A as a transition matrix to define new data Z:

.

We discovered two things. One, the matrix A is the cross covariance between Z and X:

.

I find this interesting, and I suspect that it would jump off the page at me out of either Harman or Jolliffe; that is, I suspect it is written there but it didn’t register.

Two, we discovered that we could find a matrix Ar which is the cross covariance between Zc and Xs. That is, the general relation

specializes to both centered data Xc and standardized data Xs, and gives us

.

In addition, however, we found a matrix Ar such that

,

using the new centered data Zc and the old standardized data Xs. Oh, the definition of the matrix Ar was:

,

where was a diagonal matrix whose entries were the inverse standard deviations of the centered data Xc.

That is, each row of Ac was weighted by the inverse standard deviation to define Ar. And it all worked out because the definition of Ar

was dual to the definition of standardized data, which weighted the columns of Xc:

.

I will remind you that I find this interesting, but I’m not convinced it is important. I cannot forgive Ar for the fact that it is not a matrix of eigenvectors. I assure you, I could be missing something about Ar.

Anyway, that is what we did.

What were the questions?

1. when is A a covariance matrix?
2. is v a covariance matrix?
3. can we do something similar for constant row sums?

For question 3, I have the feeling there’s something we can do, but I haven’t put the pieces together. It can’t be same as for question 1, because a matrix that weights the rows of the data matrix is not the same size as one that weights the columns (in general). Perhaps we need to look at XXT instead of XTX. But then things don’t drop out if A and Z are defined the same way, so…. This is still open.

For question 1, two partial answers are:
If , i.e. A is formed by weighting the columns of an orthonormal eigenvector matrix v by the , then A is a covariance matrix of the same form as the one from which we got v.

And the construction that gave us Ar (weighting the rows of A) works for arbitrary weights.

I’ll explain.

Suppose we have a data matrix X. Suppose we compute a “dispersion matrix”, by which I mean a multiple of :

.

(Of course, is usually 1/(N-1).)

Suppose we do an eigendecomposition of c,

;

then define a matrix

,

and finally define new data using A as the transition matrix:

.

(I’ll remind you that our observations are columns of and ; and a transition matrix applied to new components gives us old components. If you need to, look here.)

I claim that

,

i.e., that A is the “cross dispersion matrix” between X and Z.

Please note, in particular, that the data X need not be centered.

The proof goes just the way it did before.

We will need a couple of things. From

,

we get

.

And also from

we get

(because v is orthogonal and is diagonal).

Then

.

QED.

(If we had started with instead of , we would have ended up with , which is a good thing, since . Things are just a little different because the cross-dispersion in general, and the cross-covariance in particular, are not symmetric.)

This result, that A is the corresponding dispersion matrix between X and Z means that we can say something about Z without ever computing it.

(I’m not particularly fond of that, but then I haven’t been trying to deal with thousands of observations and a hundred variables. I’m used to datasets small enough to hold in my hand, as it were.)

Example

We could verify this, in particular, for the raw data that we used for the previous Basilevsky computations. We knew it was true for centered data or for standardized data, i.e. that it was true for the covariance matrix or the correlation matrix; we have seen now that it is true for an arbitrary dispersion matrix

.

It didn’t depend on 1/(N-1) or on X being centered. If we do an eigendecomposition of , then .

Well, let’s do it. Here’s the raw data that we used here:

I will take , just so that it can’t be mistaken for anything else (I need some number; let’s make sure that if anything interesting happens, we see it. OK, nothing interesting will happen, not that kind of interesting, anyway.)

We compute

Now that I’ve written it once, let me emphasize: the weird expression occurs only because I needed a number for my abstract scaling factor . In retrospect, perhaps I should not have used greek letters for both abstract and numerical.)

We get an eigendecomposition, but all I want right now is , the diagonal matrix of . OK, the eigenvalues are

and is

On second thought, I really should give you the orthonormal eigenvector matrix for reference:

The matrix A is

The new data Z is given by

Finally, the cross-dispersion between Z and X is given by

And that is the same as A, as it should be.

We continue with the answer to question 1. We found that if we do an eigendecomposition of , then .

What about the matrix Ar, found by weighting the rows of A by the inverse standard deviations?

There was nothing special about the inverse standard deviations. If we weight the columns of X by a diagonal matrix , we can find a corresponding matrix Ar.

Let’s work that out. I am going to take weights (1,2,4} in place of the inverse standard deviations. That is, I take the diagonal matrix

and scale the columns of X and the rows of A:

(Perhaps now is a good time to remind you that we can scale the columns of a matrix by post-multiplying it by a diagonal matrix of the scale factors; and we can scale the rows of a matrix by pre-multiplying it by a diagonal matrix.)

i expect that

.

Compute the RHS:

and that is, indeed, .

So, the construction of Ar using the inverse standard deviations could be done with any set of weights instead, and applied to any dispersion matrix, not just the A = Ac from centered data.

In other words, I know two kinds of transition matrices that are dispersion matrices: A between Z and X, Ar (Ad in this case) between a newer Z and the older X.

I do not have an “iff” result. I can’t guarantee that these are the only kinds of matrices that work.

Question 2. Is v a dispersion matrix?

No, at least not between its new data Z and the old X, i.e. not the way A is. (I wonder if there’s a way to make it happen. Would it ever be worthwhile if it were possible?)

This time, let’s first confirm the failure numerically.

We need new data from v instead of from A. Rather than destroy the Z’s, call the new data Y. We should have

,

i.e.

,

i.e.

.

That last equality corresponds to , so we’re good.

Here’s the new data Y from old data X using the transition matrix v:

now compute , i.e. in this case,

and that is to be compared with v – but I’ve already done the algebra, and I know they should have proportional columns, so I took the ratio instead of the difference, i.e. here is , term by term:

Those should look familiar: they are, as they should be, the eigenvalues of :

So the cross-dispersion matrix between X and Y is a column-weighted version of v, and the weights are the eigenvalues.

Time to do the algebra. This time, Y = Xv. We start with the dispersion matrix…

and that’s what we got: the cross-dispersion is the column-weighted v, and the weights are the eigenvalues .

PCA / FA Basilevsky: with data

Introduction

I am looking into Basilevsky because he did something I didn’t understand: he normalized the rows of the Ac matrix (which I denoted Ar). We discussed that, and we illustrated the computations, in the previous two posts. But we did those computations without having any data. I want to take a closer look, with data.

In contrast to As and Ac, which are eigenvector matrices, his Ar matrix is not. Nevertheless, as I said, his Ar is not without some redeeming value. In fact, all three of the A matrices have the same redeeming value.

I will show, first by direct computation and then by proof, that each of these A matrices is the cross covariance between X data and Z data.

(That doesn’t mean I want to use his Ar matrix; all it means is that I learned something new about an A matrix.)

Recall the raw data of example 9. We’re about to go get the A matrices and the new data Z for both standardized and centered data.

Yes, we’ve seen all these forms of the data before, but I want this post to be self-contained.

Now, I’m going to switch to Basilevsky’s notation, using Z for F^T. If you will, we’re mixing Harman’s model

where Z = X^T, and F^T is the new data (“scores”) WRT the transition matrix (“loadings”) A; and Jolliffe’s model

where Z is the new data (“scores”)….

We start with Harman’s model

.

eliminate Z:

,

and we know that F^T is the new data corresponding to X WRT the transition matrix A .

Now we reintroduce Z, redefining it as Jolliffe uses it, for the new data; that is, we let

and get

.

If I had to do it all over again, I would probably not have used Basilevsky’s notation. This is the notation I had for the previous two posts, but i’ll admit, it’s confusing.

But for this post, as for the previous two, that is our model.

First, let’s go get the standardized and the centered data (“X”), the corresponding A matrices, and the corresponding new data (“Z”).

Basilevsky Standardized

If you are comfortable with these calculations – which we’ve seen before for this very data – you should move right along. If, on the other hand, you need a refresher, here it is. And if you need to see the calculations explained, here

I standardize the data and call it Xs:

Get the eigenvalues of the correlation matrix:

Get the SVD (Singular Value Decomposition) of Xs:

.

Form the diagonal matrix of …

Form the weighted eigenvector matrix

Get the scores FT as the first three columns of 2u (where 2 = Sqrt[N-1], N=5). This is as close to the u matrix itself as I need to get. But instead of calling these F^T, we are calling them Z:

Check it by confirming that we have

.

(We do.)

Basilevsky centered

Recall the raw data:

I will center the data and call it Xc. Get the column means…

and subtract them from each column to get Xc:

Proceeding as before, we get the SVD (Singular Value Decomposition) and display w…

Get the eigenvalues of the covariance matrix (and that’s the only change we make after switching from standardized to centered data):

Get the diagonal matrix of of the covariance matrix…

Get the weighted eigenvector matrix

Get the scores F^T as the first three columns of 2u (where 2 = Sqrt[N-1], N=5). except that we are calling them Z, as Basilevsky does.

Check it by confirming that we have

.

(We do.)

Key: the A matrices are cross-covariances

We have As and Ac:

We have new data Zs and old Xs:

and old:

Now, what is the cross-covariance between Xs and Zs? I need either

or

,

and N-1 = 4. Note that they are the transposes of each other. We compute the former:

Not symmetric. But it’s a cross-covariance, between two different variables, not the covariance of one variable.

Now recall As:

They are the same.

For centered data:

and we recall

Again, they are the same.

Now, get Ar, by normalizing the rows of Ac:

and compute the cross-covariance between Xs and Zc:

0nce again, they are the same. We have good reason, now, to believe that

,

where Z is the new data WRT A (i.e. the scores, as most people define them).

That’s what I set out to show.

We have seen it for As, relating standardized data Xs to the scores Zs; for Ac, relating centered data Xc to the scores Zc; and finally for Ar, relating standardized data Xs to the scores Zc.

Proving it

Let’s take a look at the proofs. The question in the back of my mind is: when is a transition matrix also a cross-covariance matrix? (I will confess that I do not have a universal answer; it appears that I can generalize from the Ar matrix somewhat.)

That is, given old data X and new data Z WRT transition matrix T,

.

when can we conclude that T is the covariance between X and Z?

For either standardized data or centered data X, one argument works. We have

,

where is the diagonal matrix of and V is an orthogonal eigenvector matrix of the covariance matrix c; that is, we have the eigendecomposition

We have the decomposition of X:

.

For simplicity, i’m going to assume that A is invertible (i.e. that is invertible; i.e. that all the eigenvalues are nonzero): then

and I note that

.

Let me name those three equations: we have an eigendecomposition of the covariance matrix, the decomposition of X, and the definition of A.

Now we compute the cross-covariance. There may be simpler ways, and there are certainly other paths to the same answer, but this works for me. As I said, I’m going to assume that A is invertible, because for one thing we can cope if it is not; and for another, if it is not, we would have fewer Z variables than X variables.

(the decomposition of X)

(definition of the covariance matrix)

(the eigendecomposition of c and definition of A)

(simplify)

(definition of A)

QED.

The eigendecomposition is crucial. Without it, we are left hanging at

,

Now let us try to prove it for Ar. It’s a surprisingly simple consequence of what we just did. We need only one new equation, the relationship between standardized and centered data:

,

which says that we get standardized data from centered data by dividing each column by its standard deviation; and we recall that scaling columns can be accomplished by post-multiplication by a diagonal matrix. (After all, that’s how we get the A matrix from the V matrix.)

What we actually need is the transpose:

.

Let us compute the cross covariance of Xs and Zc.

(Xs in terms of Xc)

(by the previous result for Ac!)

(by definition of Ar)

QED.

The key, then, was that Ar was to Ac as Xs was to Xc (OK, not exactly: scale the rows of one and the columns of the other). I don’t know of any other suitable relationships, but if we found one to use, we should apply it to Xc and Ac.

Let me be more specific. We never really used the fact that was a diagonal matrix of standard deviations. If we had constructed completely different data Xd by scaling the columns of centered data Xc, using an arbitrary diagonal matrix , so that

,

and if we had defined a new A matrix by scaling the rows of Ac using the same diagonal matrix:

,

then we would find that Ad was the cross covariance of Xd and Zc.

I need to look at this some more.

PCA / FA example 9: centered and raw, 3 models

What follows is simple computation, solely to show us exactly what happens. It continues the work of the previous post, which did my default calculations for standardized data. Here I do the same calculations for centered and raw data.

Centered

The raw data is still

I will center the data and call it Xc. Get the column means…

and subtract them from each column to get Xc:

Proceeding as before, we get the SVD and display w…

Get the eigenvalues of the covariance matrix (and that’s the only change we make after switching to centered data):

I’m getting fond of just doing pie charts:

Get the diagonal matrix of square roots of those eigenvalues of the covariance matrix…

Get the weighted eigenvector matrix A

Get the scores as the first three columns of 2u (where 2 = , N=5).

Confirm that we have . (Yes.)

We compute the new data (Y = Xv = uw)

And we confirm the product – which amounts to confirming the SVD. Not a bad thing to do considering how many u’s and w’s I have floating around in my Mathematica® notebook.

Then, on second thought, I’m going to forget about computing Davis’ loadings and scores.

So we should display the two forms of loadings A, v together:

We should display the corresponding scores F^T, Y:

And we could confirm that the new data Y have mean 0 (yes) and redistributed variances equal to the eigenvalues:

(yes, the same as the eigenvalues .)

And we could confirm that F^T is standardized: the means and variances of the columns are 0 and 1 respectively. in fact F^T is uncorrelated: its covariance matrix is the identity.

What about reconstituted data? I would look at w in preference to since it’s the w’s I will be setting to zero.

As before, setting to zero will have a very small effect on the data: the reconstituted matrix would differ from the original by a square root sum of squares of 1.55 . So I’m going to set all but the first w equal to zero, just to get reconstituted data that differs nontrivially from the original.

The reconstituted data is , i.e.

and the leftover is X1 = X – X0 is

I’ll show you again that the difference between X and X0 is the sum of squares of the w’s we omitted. the squares of the entries in Xc1 are

Their sum is 27.5582 and the square root is 5.24959 . And the squares of and are 25.1575 and 2.40069, with matching sum 27.5582 .

We could compute the reconstituted new data Y, but as before, it would just be the first column of Yc.

Raw

I use the raw data and call it Xr.

Get the SVD and display w…

Neither the correlation matrix nor covariance matrix is appropriate, so we look at the singular values. But we want to honor

,

I could compute the matrix from the w matrix. Given the sizes of the matrices, I might actually compute

and take the square root. But conceptually it is cleaner to compute the eigenvalues of , which is exactly the correlation and covariance matrices for standardized and centered data. (Let me put that another way: if I did not have built-in Mathematica® commands Covariance and Correlation, I would be computing in every case.) Note, however, that for the raw data it is not remotely a covariance matrix. I get

and here’s the pie chart of them:

and then I form the diagonal matrix of square roots:

Get the weighted eigenvector matrix

Get the scores FT as the first three columns of 2u (where 2 = , N=5).

Check by confirming that for Xr, Fr, Ar. (Yes.)

we compute the new data (Y = u w):

And we confirm the product – which amounts to confirming the SVD. Not a bad thing to do considering how many u’s and w’s I have floating around.

Now let’s look at the new data Y. Here are its means…

and its variances…

The sum of the variances is 1096.36, and that is, in fact, the sum of the variances of the raw data. we have redistributed the variance.

But.

Recall the eigenvalues:

The first “but” is that the redistributed variances have nothing to do with the eigenvalues. Hey, Xc is not centered, so has nothing to do with the variance of the raw data.

The second “but” is that the redistributed variances of the centered data were

.

and the first one is larger than what I got for Yr. That’s actually good, since the first variance of Yc was supposed to be the largest I could get. The redistribution of variance of the raw data is not maximal. It’s not that far off, but it’s not maximal.

This has implications for regression analysis. That is, if I want to redistribute the variance of a bunch of variables for regression, to do it properly I have to center the data, which is not something I like doing under the circumstances. The uncentered data seems more appropriate.

The change to Y is far more significant than the change to F^T. In this case, it is still true that

,

but it is no longer appropriate to say that the F^T are uncorrelated. In fact, the means are nonzero

and the variances are not 1:

The matrix F^T is almost orthonormal, except that instead of

we have

What about reconstituted data? I would look at w in preference to :

As before, setting to zero will have a very small effect on the data: the reconstituted matrix would differ from the original by a total sum of squares of 25 . so I’m going to set all but the first w equal to zero, just to get reconstituted data that differs nontrivially from the original.

Here’s w0:

Here’s the reconstituted data:

It is still true that the raw data X and the reconstituted data X0 differ by the square root of the sum of squares of the two omitted w’s.

Finally, if we were to compute the new components of the reconstituted data we would get the first column of the Yr matrix.

PCA / FA example 9: standardized data, 3 models

Introduction

edited 16 Jan 2009: I found a place where I called F^T the loadings instead of the scores. That’s all.

I want to run thru what is admittedly a toy case, but this seems to be where I stand on the computation of PCA / FA.

Recall the raw data of example 9:

Get the mean and variance of each column. The means are

and the variances are

We see that the raw data will differ from the centered data, and that will differ from the standardized data. Let’s do the standardized data first, because that’s what we’ve been doing most recently.

Here’s what I’m going to do. For a data matrix X

1. get the SVD,
2. get the eigenvalues (in 2 cases, that’s the correlation matrix or the covariance matrix)
3. form the diagonal matrix
4. form the weighted eigenvector matrix
5. form the loadings scores
6. form the new data Y wrt v, Y = u w
7. form Davis’ loadings
8. form Davis’ scores .

I will also keep the largest “few” singular values w
and form the reconstituted data and the residual data , where:

.

From given raw data, I’m going to do all that for the raw data, the centered data, and the standardized data. Chances are only one of those is really appropriate for any given problem, but I don’t know how to decide, so I’m going to use my really cheap computing power to see it all.

That lays out everyone’s models. For Harman and Bartholomew, we’re getting

(i.e. Z = A F)

for Jolliffe and Malinowski we’re getting

(i.e. X = R C and Y = X v)

and I’ve included Davis’s R-mode scores and loadings.

This is complicated by my doing it three times, but that encourages me to check things along the way. Two of the obvious checks are

• 
• 

but we can also check that

• Y has redistributed the variance
• is uncorrelated: covariance matrix = identity

but we have to be very careful about the last two when we work with the raw data.

Finally, I will look across loadings v, A, AR and across corresponding scores Y, F^T, SR

Gentemen, start your engines.

Standardized

i standardize the data and call it Xs:

Although we saw most of these calculations in the previous post, my focus is a little different this time. Among other things, we will see just how few things need to be computed.

Get the SVD , labeled by “s”. For now, look only at w…

Get the eigenvalues of the correlation matrix:

If I simply do a pie chart, I don’t have to explicitly compute percentages:

Form the diagonal matrix of square roots …

Form the weighted eigenvector matrix A

Get F^T as the first three columns of 2u (where , N=5). This is as close to the u matrix itself as I need to get.

Check it by confirming that we have . We do.

We compute the new data Y wrt the v basis (Y = u w = Xv):

And we confirm the product – which amounts to confirming the SVD. (That is, from Y = u w, the SVD becomes .) Not a bad thing to do considering how many u’s and w’s I have floating around.

Finally, Davis defines R-mode loadings AR and scores SR:

,

so we can compute them:

OK, he would of course drop the zero columns:

In case I’ve never pointed it out before, although Davis’ scores S^R are strange, his loadings are exactly analagous to . The difference is in using w for weights instead of .
We should display the three definitions of loadings A, v, A^R (and this is when I choose to see v from the SVD computed way back at the beginning):

Since we have no zeroes, we can display the ratios (no, I wouldn’t usually do this):

i hope there were no surprises there. Each column of A or A^R is a multiple of the corresponding column of v.

We should display the corresponding scores (“new data”) F^T, Y, S^R (where F^T and Y are new components of X wrt A and v resp., but S^R is new components wrt the reciprocal basis of A.)

Since we have no zeroes, we can display the ratios:

Again, one set of scores is a column weighted multiple of any other.

We could confirm that the new data Y are centered and have redistributed variances equal to the eigenvalues: each column of Y does have mean 0, and the variances are

and those are indeed the eigenvalues

And we could confirm that F^T is standardized: each column does in fact have mean 0 and variance 1. In fact, far more is true: F^T is uncorrelated, it’s covariance matrix is the identity.

What about reconstituted data? I would look at w in preference to :

We know by now that setting to zero will have a very small effect on the data: the reconstituted matrix would differ from the original by a root sum of squares = . So I’m going to set all but the first w equal to zero, just to get reconstituted data that differs nontrivially from the original.

There are many equivalent ways to do this, as we saw last time. What come to mind are

1. take 1 column of F and one column of A.
2. set two columns of F and two columns of A to zero.
3. reduce w to a 1×1 matrix and take skinny forms of u and v.
4. set two columns of Y to zero.
5. leave sizes alone, but reset two w’s to zero.

i’ll do the last one, (5). And we didn’t actually see the second-to-last one (4), but I’ll show it to you after we get the reconstituted data. I zero-out two entries in w:

and I compute :

Now recall our initial Xs:

and look at their difference:

Incidentally, we could square each element in Xs1

and add up all those squares: the sum is 2.26886 . That, as we’ve seen before, is also the sum of squares of the omitted w’s; that is,

.

Is it plausible to ask about reconstituted Y? yeah, that’s just part of the previous calculation. We had computed ; the reconstituted Y is

Ah ha! And sure enough, the difference between Ys and Ys0 is extremely simple:

So the reconstituted Y was the first column of Ys and the other two columns set to zero. We could have just taken it and post-multiplied by vs^T to get the reconstituted X: this was variant (4) of ways to do the calculation.

OK, we’ve seen how to run thru the calculations, for the standardized data. There’s a lot redundancy in all that I’ve computed, more than I need for any analysis of my own. Still, if I’m checking someone’s else’s work, I should have routinely gotten anything they came up with.

That’s enough for now; next time, centered data and raw data.

PCA / FA Example 9: scores & loadings

I want to look at reconstituting the data. Equivalently, I want to look at setting successive singular values to zero.

This example was actually built on the previous one. Before I set the row sums to 1, I had started with

I’m going to continue with Harmon’s & Bartholomew’s model: Z = A F, Z = X^T, X is standardized, A is an eigenvector matrix weighted by the square roots of the eigenvalues of the correlation matrix of X.

I want data with one eigenvalue so large that we could sensibly retain only that one. Let me show you how I got that.

Get the SVD (Singular Value Decomposition) of t1… and look at w:

OK, keep u and v – just because they’re handy – but redefine w, and compute a new data matrix using this w:

Let :

Standardize t2 to get data X (“the data”):

Get the SVD of X, , but look at w… (we’ll see u and v later)

Get the eigenvalues of the correlation matrix and look at the percentages:

so our first eigenvalue is 81% of the sum, the second is nearly all the rest, 19%.

Get the diagonal matrix of square roots, :

Get A and the scores , except that I never want more than 3 columns. (The first 3 columns of 2 u are the components of the new data wrt the A basis.)

Check by confirming that , i.e. that we have factored the data matrix.

Let’s quickly compute the reconstituted data from 1 and 2 singular values. I could use square forms of w (w1 and w2) with u and v cut down to u1 and v1 (or u2, v2) to be conformable. To be specific,

and compute :

or take

and then :

You might note that my “1″ and “2″ indicate how many singular values I retained.

But what I was saying about leaving u and v untouched is that I can keep u and v full-size

and use the following matrices for w, all the same size as X. The full-size w is:

and the two others are:

It’s easy enough to reconstitute the data that way,

and we do get the same answers for X1 and X2.

Now, what if take just the first two columns of F^T, and the first two columns of A (i.e. the first two rows of A^T)? (We saw this last time, too.) That is, I cut F^T down to:

and I cut A^T down to

Do we get X2? Yes.

What if take just the first column of FT, and the first row of A^T? That is, I take FT to be:

and AT to be

Is their product equal to X1? Yes.

What we just saw is obvious in retrospect, but worth stating explicitly. When we throw away small singular values or eigenvalues, we do not change the scores and loadings. What we change, instead, is the number of scores and loadings used.

When we throw away a nonzero singular value or eigenvalue, we are throwing away a scores-and-loadings pair. We don’t change any of the other scores and loadings.

We have 3 pairs of “scores & loadings”, and then for each pair we compute a product. The individual pairs, and the 3 products, are not affected by our decision to “drop a factor”. Instead of adding all three products, we may choose to add only the first two products, or to keep only the first product.

It is probably customary to say that we are retaining 1 or 2 factors when we retain 1 or 2 scores-loadings pairs.

That “the scores” are selected columns of tells us that the individual scores & loadings cannot change no matter how many factors we keep.

Throwing away – choosing not to use – a nonzero scores-loadings pair does affect the reconstituted data. X2 is close to X, but not the same, and X1 is quite different.

To look at that another way, the means of X2 are zero, and the variances of X2 (2 factors retained) are:

“Not far from 1″ is an understatement.

The means of X1 are also zero, but the variances of X1 (1 factor retained) are:

You might note that the X1 data is still centered, but it is no longer standardized. That was lost when we threw away a fairly large singular value.

Next, I think I will show you what I reckon I would do, today, for PCA / FA.

PCA / FA Example 8: the pseudo-inverse

introduction

Recall Harman’s or Bartholomew’s model

Z = A F

with , X standardized, and A a -weighted eigenvector matrix, with eigenvalues from the correlation matrix.

We saw how to compute the scores in the case that A was invertible (here). If, however, any eigenvalues are zero then A will have that many columns of zeroes and will not be invertible.

What to do?

One possibility – shown in at least one of the references, and, quite honestly, one of the first things I considered – is to use a particular example of a pseudo-inverse. I must tell you up front that this is not what I would recommend, but since you will see it out there, you should see why I don’t recommend it.

(Answer: it works, it gets the same answer, but computing the pseudo-inverse explicitly is unnecessary. In fact, it’s unnecessary even if we don’t have the Singular Value Decomposition (SVD) available to us.)

Suppose we are given raw data which has, in fact, constant row sums (in this case, 1):

A couple of those fractions do not display well on my screen, so let me convert to decimals:

PCA a la Harman or Bartholomew et al.

Now, we do a PCA using Harman’s or Bartholomew’s model – but my way, not theirs. Standardize the input:

The X matrix is “the data”.

Get the eigenvalues of the correlation matrix:

Yes, one of them is identically zero. We learned here that this happens when we start with constant nonzero row sums and compute the correlation matrix.

Get the diagonal matrix of square roots:

Get the SVD of X… and look at w:

Get the -weighted matrix A, as . (The only reason I’m using “1″ is so that in subsequent algebra and computations I can use A, v, etc. consistently for cut-down matrices.)

Compute as I would usually, as columns of , in this case 2 u1, and only the first 3 columns:

Let’s check it by computing :

and that is, indeed, X:

We’re done.

But what if the SVD is not available to us, and we don’t have u?

Let us recall the derivation in the case that A was invertible (here it was called A1; this is old stuff, so try not to worry about all the missing “1″s), we could also have computed by projecting X onto the reciprocal basis :

,

but we can also derive that formula without ever knowing about the reciprocal basis:

.

In that case, we could also have used to simplify things further, and get

(because v is orthogonal)

and finally

.

That was fine (and it’s what Bartholomew had us do), but our A is not invertible.

the pseudo-inverse

So let’s work with the cut-down matrices. To be specific, we start with the two nonzero columns of A1:

This cut-down A has, in a sense, eliminated the problem: although the new A can’t be inverted because it’s not square, it is of full rank. Even better, it’s A1 with one column removed: in a very real sense, it’s the same linear operator.

Any time I have a formula which involves a matrix inverse, and I want to generalize it to a case where I have a rectangular matrix – hence, no inverse – I investigate whether a pseudo-inverse will work.

In this case, A is of rank 2, so the two square matrices and are of rank 2. But A is 5×2, so is 2×2 and is 5×5, so is invertible, even though A1 wasn’t. There’s no point to considering , because it’s not invertible.

So, for the rectangular A – because what I dropped was a column of zeroes – we still have

Z = A F

(That is crucial. We didn’t change the data , because we didn’t lose any singular values.)

Premultiply by :

.

now premultiply by , which we know exists:

.

It may not be pretty, but it will work.

I should remind us that this is exactly what’s goes on in regression (OLS) here: the “normal equations” for the regression model

are

.

Let me also point out that if A is invertible, then the inverse expands as

,

so from

we get

i.e.

and finally

.

That is, the equation written with the pseudo-inverse

contains our earlier special case

when A is invertible. That derivation, incidentally, reminds me that although I think of as the pseudo-inverse, it’s really the product which is the pseudo-inverse, since it’s the product that collapsed to .

It gets better, in this case. The pseudo-inverse is far prettier than it looks, but let’s see it before we show it.

Here is :

Diagonal?

Diagonal!

Even better, it’s , where is the cut-down :

That is, it’s the diagonal matrix of instead of . Not only is it diagonal, but its elements are almost already computed. And in fact the inverse is almost immediate:

.

Then a quick check shows that just as

,

we have

,

hence

(We also need a cut-down version v of v1 (its first two columns) as well as the cut down of (2×2) which we got earlier.)

Then is

(That is, whether we start with A1 and cut it down to A, or compute the cut-down A from v and , we get the same thing.)

Putting it all together, from

we get, as before,

.

That is, if A is invertible, we can compute as any one of

• 
• 
• (the first “few” columns).

If A is not invertible but the cut-down A is of full rank (so v and are also cut down as above), we can compute as any of

• 
• and take the transpose.
• (the first “few” columns).

The only difference between A not invertible and A invertible is in place of .

(Yes, one of those uses F, the other .)

But as I said at the beginning, I didn’t go thru all that in order to encourage you to compute the pseudo-inverse , but unfortunately you may very well see it in texts.

When you do, realize that it’s true but utterly unnecessary for computing F.

(There’s no reason to tell you which of the references show it.)

I would always compute as . If the SVD makes you uncomfortable, or isn’t available (I weep for you), then compute it as , using V from the eigendecomposition.

Even though the pseudo-inverse isn’t all that useful in this case, it’s a good thing to know about in general.

Oh, and don’t even dream that it’s unnecessary for regression. For PCA / FA, it’s

which makes it unnecessary to actually compute the pseudo-inverse . For regression, isn’t at all likely to be diagonal.

PCA / FA tricky preprocessing

Introduction

I have stumbled across a tricky point in the preprocessing of data. The most relevant post is probably

this of April 7. Rather than lecture, let me ask and answer some questions. The fundamental question is:
Can I inadvertently reduce the rank (the dimensionality) of the data matrix?
The answer is yes.

Suppose we have a data matrix of full rank. Suppose the matrix is a typical one, in having more rows than columns. That the matrix be of full rank is equivalent to: the columns are linearly independent.

There is a serious asymmetry here: not all of the rows are linearly independent.

We might as well imagine some specific numbers. Suppose we have 3 columns and 5 rows. Since we assume it is of full rank, the rank of the matrix is 3.

The term “centering” means to “make zero-mean.” A centered column is one with zero mean = zero sum; a centered row is a row whose mean (hence sum) is zero..

What happens to the rank of the data matrix if we only row-center the data?
The rank decreases by 1.

Visualize the matrix as 3 column-vectors c1, c2, c3 of length 5. That the row sums (= row means) are zero says that c1 + c2 + c3 = 0, so the three columns are linearly dependent and the matrix is of rank 2.

What happens to the rank of the data matrix if we only column-center the data?
Nothing.

Visualize the matrix as 5 row-vectors r1, r2, r3, r4, r5 of length 3. We have more than 3 vectors of length 3, so they span a 3D space, but no larger. Adding the constraint r1 + r2 + r3 + r4 + r5 = 0 imposes, in general, no additional restriction. (If they had somehow spanned a 5D space, it would, but they can’t.)

What happens to the rank of the data matrix if we only set the row sums to 1 (or 100), i.e. to a nonzero constant value?
Nothing.

All we’ve done is assert that r1 + r2 + r3 = 1, a constant vector. These 4 vectors are linearly dependent, but the original 3 are not.

But. It has another potential effect, probably undesirable.

If we were to use the resulting data matrix for a regression, we would find ourselves in trouble. We almost always want to adjoin a column of 1s to the data matrix, in order to include a constant term in the regression equation. We just got through saying that the combination r1, r2, r3, 1 is linearly dependent: the regression will explode with a singular matrix. Once we have set the row sums, we would have to drop a column of data in order to include a constant term. This is the same problem we face in regression when we define instrumental (dummy) variables.

What happens to the rank of the data matrix if we only set the column sums to 1 or some other nonzero constant value? (I’ve never seen it done, but why not, at least in principle?)
Nothing. (It even has no effect on regression.)

First Summary

So, there were 4 cases: whether we set row or column sum to a constant value, and whether that constant was zero or not. Of those 4 cases, only one (row sum = 0) reduces the rank of the data matrix.

What I was unclear about – or downright wrong about – on April 7 was in saying that we would lose rank by setting the row sums to 1 or 100 etc. Not so; none of the other 3 cases reduce the rank of the matrix.

Is that it?
No. All those questions were about one single operation. What if we combine operations?

Can’t we center both the rows and the columns?
Yes. We saw
here on April 13 that we could doubly-center the data.

Does that mean we lose rank?
Yes, because we are centering the rows. Centering the columns is irrelevant.

Can we set the row sums to 1 and then center the columns?
Yes.

But something else will happen. Centering the columns is no longer irrelevant.

Didn’t we see that if the row sums were constant, and then we centered the data, the row sums went to zero?
Yes. It may seem weird, but it’s true.

That’s an easy way to get doubly-centered data using Mathematica®. Scale the rows so they each sum to 1, then center the columns. Voila’, the row sums will change to zero.

Any time you start messing with row and column sums or means, remember the grand mean, too: the average of all the data. It is also both the average row mean, and the average column mean. If the columns all have zero mean, then the grand mean is zero. Then the mean of the rows (edit: i.e. of the row means) must be zero. The challenge is to show that the rows still have a common sum (edit: equivalently, a common mean); then that common sum mean must be the mean of the rows, i.e. zero.

If we set the row sums, say to 1, and then column-center the data, do we lose rank?
Yes.

Column-centering (or standardizing) the data after setting the row sums will reset the row sums to zero, and that loses rank. If it is appropriate, however, that the data be column-centered, so be it. I do not believe that reduction of rank is sufficient grounds for not centering the data.

Second summary

The apparently innocuous act of column-centering can lose rank, if we apply it to nonzero-constant-row-sum data. We can do it, and sometimes we should do it, but we should expect to lose rank when we do. But we’re losing rank because we first set the row sums, even though to a nonzero constant.

I had previously said I would favor doubly-centered data. Now I am inclined not to do so without cause, simply because it reduces the rank of the matrix.

Can I inadvertently reduce the rank (the dimensionality) of the data matrix?
Hell yes. In particular, we might find that we effectively centered the data without realizing it.

Suppose we set the row sums of our data to 1, say, and then we decide that rather than do an SVD of the new matrix, we will do an eigendecomposition of the covariance (or correlation) matrix.

Oops. Computing either the covariance or correlation matrix is tantamount to column-centering the data.

here on March 10.

Goodbye full rank.

Third summary

Suppose we have data for which a constant nonzero row sum is appropriate. Suppose we choose to do an eigendecomposition of the covariance (or correlation) matrix.

Then we lose rank. Constructing either of those smaller matrices implicitly centers the columns, and that in turn implicitly changes the row sums of the data to zero. Even if we never explicitly compute the transformed data.

The clearest path in such a case might well be as follows. First set the row sums to 1 or whatever nonzero. Then column-center or standardize the data. We will see the rank drop by 1 at that point. Then do an SVD. The only advantage to doing it this way, rather than via the covariance or correlation matrices, is that we see when we lose full rank.

If we live in a universe where data is supposed to be centered, or where the covariance or correlation matrix is appropriate, then we will lose rank if, in addition, the data had constant row sums.

Final summary

0. We could set row or column sums to zero or to a nonzero constant.
1. Setting row sums to zero is the only one of the four cases that will reduce rank.
2. Column-centering applied to data with constant row sums leaves us still with constant row sums but the sums have been changed to zero. This gives us (1) and loss of rank.
3. An accidental way to achieve (2) and thus (1) is to analyze the covariance (or correlation) matrix of constant-row-sum data.

It’s not a problem if you know what happened and why.