Example: Is it a transition matrix? Part 2

We had three matrices from Jolliffe, P, V, and Q. They were allegedly a set of principal components P, a varimax rotation V of P, and a quartimin “oblique rotation” Q.

I’ll remind you that when they say “oblique rotation” they mean a general change-of-basis. A rotation preserves an orthonormal basis; a rotation cannot transform an orthonormal basis to a non-orthonormal basis, and that’s what they mean — a transformation from an orthonormal basis to a non-orthonormal basis, or possibly a transformation from a merely orthogonal basis to a non-orthogonal one. In either case, the transformation cannot be a rotation.

(It isn’t that complicated! If you change the lengths of basis vectors, it isn’t a rotation; if you change the angles between the basis vectors, it isn’t a rotation.)

Anyway, we showed in Part 1 that V and Q spanned the same 4D subspace of .

Now, what about V and P? Let me recall them:

As before, I suppose there is a transition matrix T relating them; each column of can be computed by applying a transition matrix T to the corresponding column of :

V’ = T P’

i.e.

V = P T’.

Then I imagine briefly that I can “solve” for T’…

,

but then quickly replace the non-existent inverse of Q by the appropriate pseudo-inverse

.

I compute T’, and then display its transpose T:

and then I test whether

V = P T’ (?).

The largest difference (element by element) between the matrices V and P T’ is 0.170627, which is rather larger (by a factor of 23) than the differences we saw in Part 1. Here’s the full blown matrix — where I chose to round off first, then subtract (i.e. what we’d compute if we were subtracting published, already rounded-off, matrices):

I am inclined to think there is a probem. But we don’t know much at all about it.

Let’s try one of our alternatives. We construct a projection operator onto the subspace spanned by P.

Get the SVD (Singular Value Decomposition) of P…

.

Then, since u is 10×10 and V is of rank 4, the leftmost 4 columns (u1) of u….

are an orthonormal basis for the range of P (i.e. the column space of P). We construct a projection operator, as we have before, by

.

We check it first by confirming that R R = R, i.e. it is idempotent, hence a projection operator. We check it further by applying it to the column vectors in P: we expect that R P = P; i.e. it actually projects onto the column space of P (the space spanned by the columns of P).

Now, having checked that R is a projection operator onto P, we apply it to V, and compute V – R V:

Let me round that to the nearest .01…

I should put that array of differences into perspective. Recall V:

This reinforces my belief that we have a problem. Some of the differences we found are larger than some of the values of V itself. I can see that V is not in the column space of P, and while it’s not extremely far away, it’s not close either. (This seems more definitive than “the hypothetical transition matrix T doesn’t quite work”.)

As before, however, we are tip-toeing around the real question: Is every column vector in V a linear combination of the columns of P?

Easy enough. Do a regression. Actually 4 of them, one for each column of V.

Let me name the individual columns of the P matrix, and let me set the first column of V (v1) as “y”, and the columns of V as independent variables.

(Some smaller pictures follow, in which you can see that the final option is “IncludeConstantBasis -> False”.) The results are:

Not bad at all. In fact, rather close to perfect.

What was the first row of T?

{0.905145, -0.386677, 0.0882059, -0.141151}.

Yes, those are the (“Estimates”) for that regression; first row of T is first column of T’ is regression coefficients for the first column of V as the dependent variable.

Let me emphasize another relationship. For the first regression, yhat should be the projection of v1 (the 1st column of V) onto the column space of P. But we have that projection directly as R v1:

R v1 = {0.477784,0.476228,0.355644,0.257779,0.487767,
0.058964,0.209419,0.0529866,0.123281,0.19278}

… and we have the yhat from the regression:

yhat = {0.477784,0.476228,0.355644,0.257779,0.487767,
0.058964,0.209419,0.0529866,0.123281,0.19278}

Good, they are the same. So if regression works better for you than the projection operator, then use regression. Personally, I prefer the projection operator because I can apply it to all of V at once. (And because sometimes I have a use for the orthonormal basis u1.)

Largely to display the R^2, let me run the other three regressions. (Oh, I also reduced the magnification in my Mathematica notebook, so the pictures are smaller and the right side of the regression command isn’t cut off, at least on my screen.)

Here’s v2 as the dependent variable:

Here’s v3 as the dependent variable:

Here’s v4 as the dependent variable:

We see that the is over 99% for v2, falls to 95% for v3, and is 97.6% for v4. Those last two are just not large enough, not for this question.

I’m certain that V and P do not span the same subspace, but that’s about all I know.

We could, however, have made the dual choice, projecting P onto V instead of V onto P. Are they the same?

No.

We already know — from the T matrix, from the projection operator, or from the regressions — that the P and V matrices do not represent the same data in two different coordinate systems. But there might be — and in this case, there is! — more information available by making the other choice.

Let’s get the projection onto V instead of onto P. We need the SVD of V…

.

We need the 4 leftmost columns (u1) of u…

We construct the projection operator…

Now we apply it to P, and compute the difference P – R P:

So what?

Round it off:

Ah ha! It’s time for Robin to squeak out, “Holy moly, Batman!”

The fourth column of P does not lie in the subspace spanned by the 4 columns of V, while the first three do. That is the problem. I don’t why or where it happened, but I know exactly what’s wrong: the 4th column of P is the only problem.

Somebody blew it somewhere. (OK, since I don’t know how a varimax rotation is computed, I have to entertain the possibility that V is not supposed to be P under a change of basis. But I can’t believe that, because V and Q have a common column space, and 3 of the 4 columns of P lie in that common column space. If that’s an accident, what does a deliberate act look like?)

Is the mistake in the orignal paper? Were P and V computed from the same correlation matrix but by independent algorithms? Or was V computed incorectly from P (and then Q computed from V)?

Or did someone copy P incorrectly? (P should be primary, but then how the heck does the 4th column of P fail to be in the column space of V? I would sooner have expected an error in V, given P.)

Or did they rotate the first 3 PCs, and then apply that transformation to the 4th column? (No. I’ve checked. You can too.)

Let me emphasize that we would get the same result from four regressions each with one column of P as the dependent variable and all the columns of V as the independent variables. The first three regressions would be extremely good fits, but the fourth would not be. We would, again, conclude that the fourth column of P was not exactly — or even to a fair approximation — a linear combination of the columns of V.

Summary: theory

Given two matrices, say P and V, which are alleged to be the same data in different coordinate systems, we may check that assertion very quickly (if both are of full rank):

Define

.

If we find it true that

V’ = T P’,

then T is a transition matrix between two bases, and P and V are indeed the same data in different coordinate systems. We’re done.

If, however, the equality fails, ,

then T is not a transition matrix, and P and V are not the same data in different coordinate systems. Unless we want more information, we’re done.

Finally, we have seen that, numerically, we may have T “close to” a transition matrix. (In Part 1, V and Q are very nearly the same data; in the color posts, the xyx bar and rgb bar tables are very nearly the same data.)

Summary: practice

Computing T tells us if something is wrong, but not much more. It’s quick, easy, and gives a simple answer: “good” or “no good”.

Computing both projection operators P onto V and V onto P (not just one of them) would likewise tell us whether or not something is wrong, and, in addition, the projection operators might isolate the problem.

We could use regression instead of computing explicit projection operators.

Summary: this example

We found that V and Q are the same data in two different coordinate systems, but that P is not. In addition, we actually know that the problem is limited to the 4th column of P: it does not lie in the space spanned by the columns of V (i.e. it does not lie in the column space of V, which is also the column space of Q).

Question:
Can someone tell me if these P, V and Q match the original paper by Yule et al.?

And there’s a question you might have.
How many examples did I have to search through to find such a good (because it was such a bad) example?

One. This was the very first case I checked, outside of the 1931 CIE tables.

That’s scary.

Example: Is it a transition matrix? Part 1

This example comes from PCA / FA (principal component analysis, factor analysis), namely from Jolliffe (see the bibliography). But it illustrates some very nice linear algebra.

More precisely, the source of this example is:
Yule, W., Berger, M., Butler, S., Newham, V. and Tizard, J. (1969). The WPPSL: An empirical evaluation with a British sample. Brit. J. Educ. Psychol., 39, 1-13.

I have not been able to find the original paper. There is a problem here, and I do not know whether the problem lies in the original paper or in Jolliffe’s version of it. If anyone out there can let me know, I’d be grateful. (I will present 3 matrices, taken from Jolliffe; my question is, does the original paper contain the same 3 matrices?)

Like the previous post on this topic, this one is self-contained. In fact, it has almost nothing to do with PCA, and everything to do with finding — or failing to find! — a transition matrix relating two matrices.

On p. 163, Jolliffe provides a matrix of “principal components”:

On the same page, he provides two matrices of the “rotated factor loadings”. One is “Varimax”…

… and the other is “direct quartimin”:

I presume — not knowing how to comute a varimax rotation! — that each of these matrices represents the same data, i.e. the same 4D subspace of R^{10} (since we have 4 column vectors of length 10 in each matrix). It would seem silly if a varimax rotation of PCs was not supposed to represent rotated PCs within the same subspace.

That presumption is wrong: they do not all represent the same subspaces.

How do I know this?

Because given two matrices of the same shape, I know how to find the transition matrix if it exists. We did this here.

Let me elaborate on what I’m doing. Imagine that our matrices had 2 columns and 3 rows. Those 3D column vectors lie in a 2D subspace — a plane — in . By doing 2D transformations within that plane, I might find a particularly nice basis for representing that data.

If, instead, I apply a 3D transformation, I would be moving that data into a different plane. What I am doing here is looking at 4D transformations within the 4D subspace. Everything is predicated on my belief that the PCs in P, and the varimax rotation of them in V, and the direct quartimin data in Q are all supposed to be in the same subspace.

Not knowing how they were supposed to be computed, I could be seriously wrong here. OTOH, since we will find that V and Q are the same subspace, and the first 3 of the 4 columns of P also lie in that common subspace, I really believe that the 4th column of P should too. That it does not strikes me as an error.

Let me make one simplification up front. Are the matrices P, V, Q of full rank (4)?

Yes. The Mathematica command MatrixRank says that each of those matrices is of rank 4.

Let us consider V and Q. If the observations in each are related by a transition matrix T, and if we arbitrarily take V to be the “old” data, then each column of V’ (the transpose of V) is found by applying the transition matrix to the corresponding column of Q’.

V’ = T Q’

i.e.

V = Q T’.

Then, as I showed before, I imagine briefly that I can “solve” for T’…

,

but quickly replace the non-existent inverse of Q by the appropriate pseudo-inverse

.

(It’s worth repeating that the inappropriate pseudo-inverse would use QQ’, but that can’t work because QQ’ is 10×10 but of rank 4 at most, hence not invertible.)

If Q’Q is invertible, that recipe always gives me a matrix T, but T need not be a transition matrix. The resulting T is a transition matrix if and only if we actually have

V = Q T’,

(given that V and Q are of the same full rank).

Do we have V = Q T’ ?

By direct computation I get

and then I look at V – Q T’…. You know, it’s easy enough to put them side by side.

Even the difference of the rounded-off matrices is satisfying small. That is,
Round[V,.01]-Round[Q.TT,.01]//Chop//MatrixForm is

Judging that the difference is negligible, I conclude that the matrices V and Q do span the same subspace of ; they are the same 4D subspace with two different coordinate systems (two sets of basis vectors).

While we’re here, let’s do this another way. After all, the fundamental question is: can the columns of V be written as linear combinations of the columns of Q? Alternatively, can each column of Q be written as a linear combination of the columns of V?

Bear in mind that we have already established that Q and V are of the same full rank. In that case, the alternatives are equivalent.

We can answer the fundamental question by trying to regress each column of V on Q. Just remember to tell the regression to drop the constant!

To be specific, here’s the first column of V as a linear combination of the 4 columns of Q. Here’s the parameter table and — not the adjusted : I really want to know how good this fit is, never mind how many variables it has.

Here’s a snapshot of the Mathematica commands (version 7). (Oh, darn. The last option in the LinearModelFit command, which replaced “Regress” is “IncludeConstant->False”. And I print the names just to confirm that none of those symbols has a left-over numerical value!)

(Let me just remind you that the is a measure of a how close the fit is. When we add another independent variable, however, the cannot decrease; it is not a good measure of the value — i.e. the worth — of that latest independent variable. The adjusted can decrease; in fact, in my experience, if the t-statistic of that latest variable is greater than 1 (in magnitude), then the adjusted goes up. As I said, in this case, I really want the rather than adjusted because I know exactly how many independent variables I am supposed to use. That said, there isn’t a whole lot of difference between the two in this case.)

That’s a really good fit. And it should be. And that column of , i.e. “Estimate”? It’s the first column of T’, i.e. the first row of T:

{0.9273, 0.0921888, 0.15161, 0.101824}

We could run 3 more regressions, but they’d only confirm what we believe: T’ could, in fact, be computed as regression coefficients. If the two subspaces coincide, then T is also a transition matrix between two bases in one and the same space.

There is yet another way to look at that. We know that in regression, yhat is the projection of y onto the subspace spanned by the columns of X. If y itself is in that subspace, then yhat = y, and the relatonship is given by a transition matrix.

So we could do this a third way. We know how to compute a projection operator onto the subspace spanned by the columns of X. In this case, X is Q, so we find its SVD (Singular Value Decomposition)…

Then the 4 leftmost columns (u1) of u are an orthonormal basis for the range of Q….

That is, u1′ u1 = I, if the column space of X is of full rank. Check it:

Bur the reversed matrix product

R = u1 u1′

gives us a projection operator R onto the range of Q:

We can check it. (I won’t print it: it’s 10×10.) One, R is a projection operator if and only if it is idempotent, i.e. R R = R. And by computation, I see that is true.

Two, applied to each column of Q, it reproduces that column of Q; i.e. applied to Q as a whole, it reproduces Q, i.e. R Q = Q. Again, computation shows that is true, too.

Having checked that R is a projection onto the column space of Q, let’s use it: apply it to V. If V lies in the subspace spanned by the columns of Q, then R V will reproduce V….

I find that the largest difference between the elements of R V and V is 0.00720726 .

That’s as close as the other two calculations. Not exactly zero, but then, our regression didn’t fit exactly: like the RGB and XYZ tables in the earlier post, V and Q are almost but not exactly related by a transition matrix.

So we have seen three ways of looking at this. All three tell the same story: the column spaces of V and Q are the same.

Let me say a little more about two of them, the projection operator and the explicit regression.

The projection operator R came from the SVD. The regression can also be written using a projection operator, namely the hat matrix H. I keep getting mixed up about the relationship between T and R and H, so let me write it out.

The normal equations for

are

.

(I showed you a quick and dirty way to recover that equation in the earlier post.)

Then

and so we may write

yhat = H y

with the hat matrix H defined as

.

For our regression, we took y to be the first column of V, and X to be the four colums of Q, so our hat matrix would be

.

It is a projection operator onto the range of Q. But that’s exactly what R was. Go ahead and check: you will find that

H = R.

This means that we could compute R as the hat matrix, instead of using the SVD. Personally, I know how to use the SVD, whereas I would have to work out the definition of the hat matrix every time. But that’s a preference; you may choose to compute the hat matrix instead. (And that is how Malinowski was doing it!)

Going a little further, T’ was

so we have

Q T’ = H V.

We’ll look at the relationship between P and V in Part 2.

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. Example 4 ! Davis, and almost everyone else

I would like to revisit the work we did in Davis (example 4). For one thing, I did a lot of calculations with that example, and despite the compare-and-contrast posts towards the end, I fear it may be difficult to sort out what I finally came to.

In addition, my notation has settled down a bit since then, and I would like to recast the work using my current notation.

The original (“raw”) data for example 4 was (p. 502, and columns are variables):

It would be prudent to compute the means…

and variances of the data…

It would also be prudent to calculate the row sums.

They are all the same (we have “constant row sums”).

Now, Davis chooses to analyze centered data. Let’s get it, by subtracting its mean from each column.

That is what he has on page 503. I remind us that we expect (HERE) that the row sums of the centered data are now zero. That in turn means that we have lost rank: the centered data should have two nonzero eigenvalues or singular values, instead of 3.

Next, Davis chooses to analyze instead of the covariance matrix . We recall that these two have the same eigenvectors, but eigenvalues differ by a factor of N-1. In fact, Davis also uses the singular value decomposition (SVD); that is a good reason for doing an eigendecomposition of .

Here is the SVD of X, . (Always remember that my u and v correspond to Davis’ V and U respectively, if you’re looking at both of us.)

Sure enough, there are only two nonzero singular values in the w matrix.

Let us now compute Davis’ loadings and scores. Using my notation for the SVD, he has

, but use the definition for computation: one less matrix product.)

(I should remark that we see that we could have done all of those calculations using the data X, the singular values w, and the matrix v: we do not need u.)

and we get

Davis does not have any of these extra columns of zeros, but all of our nonzero numbers agree (to within a sign for each column). At this point, I think the most significant thing to say is that these scores and loadings do not reproduce the data; in the sense that, the matrix product of corresponding scores S and loadings A is not X.

Let us see what Jolliffe would have done.

He writes his model as

Z = X A,

with Z as the new data, and A as an orthogonal eigenvector matrix; but I would write that as

Y = X v,

since I reserve A for the weighted eigenvector matrix.

(I have to point out that, as we saw in previous posts, .)

Now, Jolliffe would usually do an eigendecomposition of the covariance matrix; but we have enough experience to understand that we can do an eigendecomposition of . Here it is:

Here are its eigenvalues:

Here is an orthogonal eigenvector matrix (now I use my standard notation, V for an orthogonal matrix computed from an eigendecomposition, versus v from the SVD):

We can check that . (It does.)

We know that V should differ from v (from the SVD) by at most the signs of columns. Recall v from the SVD:

In fact, they agree completely.

Jolliffe would project the data X onto the eigenvector matrix V to get new data Y = XV:

Now we should recall , to confirm the equivalence:

We recall that Y = X v = X V is equivalent to the change of basis formula

x = v y,

where x is an old observation and y is the corresponding new one.

We learned – from Harman! – that the new data Y has maximally redistributed variance. It is crucial that we used an orthogonal eigenvector matrix v or V. Let’s confirm that by computing the variances of the columns of Y…

and their sum is 32.

Let’s recall the variances of the columns of X…

and their sum is also 32.

Now, what about Harman? He writes his model as

Z = A F,

where

,

(i.e. his variables are in rows, not in columns) and A is a -weighted eigenector matrix. And I generally favor his notation for Z, but I will follow Jolliffe for a few seconds later on.

Let us compute. From the eigenvalues of c,

we construct a diagonal matrix of their square roots:

We should expect that these are the nonzero nonsingular values w. Sure enough:

Note, however, that w is not square, hence is not equal to the diagonal matrix . This will lead to an irrelevant difference between A and , etc. In fact, I often ignore the distinction between and w conceptually; in practice, of course, whatever matrices we’re using have to be “the right size” for whatever matrix products we compute (conformable).

We should recall that the nonzero eigenvalues are the nonzero values of both and :

The nonzero eigenvalues were…

whereas

and

Resuming our thread, we compute the weighted eigenvector matrix A

and that should be Davis’ ; and it is, except for an extra – to my mind, irrelevant – column of zeros.

Viewing A as a transition matrix, we know that Harman’s scores are the new data . We also know (HERE) that we may compute them as

,

that is, the “first few” columns of u, not all of u.

Ah, life gets interesting because A is of rank 2 instead of 3. The new data should be the first two columns of u, but if we want to see (Harman’s) Z = A F, then we must either cut down A or keep the first 3 columns of u, to have conformable matrices. I will cut down A.

We check that by computing A F and comparing it to ; we expect Z = AF. And we have it; here’s Z:

And here’s A F (using the two nonzero columns of A!):

I remind us that this new data does not have redistributed variance; in fact, this new data has constant variance:

There are a couple of (intertwined) ways to understand why the variances are not 1. Fundamentally, we did an eigendecomposition of rather than of the covariance matrix . Alternatively, our new data was given by columns of u rather than by columns of . (And that’s a consequence of the “fundamental” reason, but it may make sense by itself.)

That constant variance is in marked contrast to the maximally redistributed variances of the new data Y constructed from the orthogonal matrix v instead of from A.

It is my strong preference to compute the new data from the first few columns of u. To make sure I have the scaling correct, I then check the product A F.

Any number of people, including Davis himself in a subsequent example, show us how to compute the new data using a reciprocal basis instead; that is, one which weights the v matrix by instead of . (Whether or not an author uses the phrase reciprocal basis or dual basis, if he divides each column of v by , he is computing a reciprocal basis.)

Life is made interesting, but not difficult, by the presence of a zero eigenvalue. Frankly, the simplest way to deal with a zero eigenvalue – for finding a reciprocal basis – it is to set the zero eigenvalue to 1. Instead of the diagonal matrix

I use

and instead of , I compute the reciprocal basis (i.e. I divide each column by instead of multiplying by ).

Then the new data with respect to A can be found by projecting X onto B:

That had better be … and it is:

(Okay, that was the simplest way to get a reciprocal basis when we have a zero eigenvalue; but the simplest way to get the new data is to avoid the reciprocal basis entirely by taking the first few columns of u instead – assuming that you have the SVD available to you.)

For this example, Davis did not compute the new data with respect to A (i.e. wrt his ). As I have said before, I have seen him do it in another example; but his scores are not the new data , which is, as far I have seen in textbooks, what other people use for the scores corresponding to loadings A.

Moving on, we learned from Basilevsky (HERE) that the weighted eigenvector matrix was also the dispersion matrix between the new data and the old data X. I say “dispersion matrix” because A itself was not computed from a covariance matrix but from a more general “dispersion” matrix .

Let’s confirm that. For just a moment I want to denote the new data by (Jolliffe’s) Z instead of , so as not to get confused by transposes. A quick computation (on the side) shows me that I want to compare A with , i.e. that I want to compute

(finally reverting to my original notation). Fine, I compute that product, and I get

and I compare that to A:

Good. Our matrix A is, as we expect, a cross-dispersion matrix between old and new data. (Perhaps I should be less meticulous, and refer to A as a cross-covariance matrix; that term may be more familiar and hence clearer in sense, even though it is not the covariance matrix.)

Now let me show you something new. We had also discovered from Basilevsky (HERE) that the orthogonal eigenvector matrix v was not a dispersion matrix D between its old data X and its new data Y. The v matrix does not have this property of the A matrix. We had seen that v and D differed by eigenvalue column-weights:

.

But… we’ve just seen that weight here. Okay, it wasn’t , but it was either or . Well, it must have been , i.e. the one that multiplied v:

.

Of course, I should check that computationally. Here is :

and here is our :

Hang on. We need to explicitly see , the cross-dispersion matrix analagous to : we compute it…

I don’t know if it’s useful, but there it is anyway. I can hear my favorite Vulcan whisper, “Fascinating.”

Don’t misunderstand me: I am sure that the Q-mode scores and loadings are interesting in their own right. (I don’t understand the goal of a Q-mode analysis, but that’s another question.)

But I do find it… fascinating… that Davis’ Q-mode scores are the cross-dispersion matrix between the old data X and the new data Y (in addition to our earlier finding that his Q-mode loadings are the new data Y with respect to the transition matrix v).

To put that another way…. The orthogonal eigenvector matrix v could have been computed from an eigendecomposition of the dispersion matrix , and defines new data Y; and is the cross-dispersion between old data X and new data Y.

I had asked a question: is the orthogonal matrix v the dispersion matrix D between the old data X and the new data Y? The answer was, “No” (HERE).

Now we know the rest of the story: not v, but the Q-mode scores are that dispersion matrix D.

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 Basilevsky: what he did.

Introduction

The previous post discussed something interesting that Basilevsky did. (Search the bibliography, he’s there.) I can’t say I like it — because it leads to a basis which is non-orthogonal, and whose vectors are not eigenvectors of either the correlation or covariance matrices.

But I had to understand it.

I don’t know how widespread it is nowadays, but even if Basilevsky is the only author who does it, it’s another example of the lack of standardization (no pun intend, I swear) in PCA / FA. This branch of applied statistics is like the mythical Wild West: everybody’s got a gun and there are bullets flying all over the place. Law and order have not arrived yet.

OTOH, it’s nice to find something different in just about every book I open.

Let me set the stage again. What we have is the following 3 models.

.

where X is data (variables in columns, observations in rows, and there are N rows), Z is principal components (PCs), and A is constructed from an eigenvector matrix V using the diagonal matrix of .

Now is a good time, before we need to be reminded of it, to emphasize that all three models are of the same form:

.

In the first model,

,

Xc is centered data, and

,

where Vc is an orthogonal eigenvector matrix of the covariance matrix :

.

In the third model,

,

Xs is standardized data, and

,

where Vs is an orthogonal eigenvector matrix of the correlation matrix :

.

In the second model,

we have the standardized data Xs as in the third model, but the PCs Zc from the first model. The transition matrix Ar may be found by normalizing the rows of Ac:

,

where is a diagonal matrix of the standard deviations () of Xc.

The example

Basilevsky gives us a covariance matrix c. (So we have no data, for this example.)

We will want both the correlation matrix r and the diagonal matrix , so let’s get them. is to be a diagonal matrix of the square roots of the variances, hence of the square roots of the diagonal of c. (i should call it s, because it’s made up of sample standard deviations, but I don’t want “s” as well as “As”, “Vs”, etc.) Here it is:

Once we have , the fastest way for me to get the correlation matrix is

:

From the covariance matrix:

So let’s do the eigendecompositions. First from the covariance matrix, which will give us the first model. We get eigenvalues…

an eigenvector matrix…

These are not quite Basilevsky’s numbers. The 4th eigenvector is off, by as much as .0014 . ah, the 2nd is off too, in two terms. Not great, but not a big deal. (Actually, I want to check his numbers, but not today. I will be using my numbers for eigenvectors and eigenvalues.)

(I suspect that most packages which do PCA / FA use iterative schemes to find one eigenvector after another; I’m getting used to seeing more and more errors as we move to the right in the eigenvector matrices. I think they’re suffering from accumulated round-off errors. Of course, I generally take it for granted that my answers are the right ones, but see below.)

Now I’m going change the signs to agree with Basilevsky:

We will discover that Basilevsky should have changed the sign of the second column, too, but he didn’t so I won’t either.

Let’s keep going. We get the diagonal matrix of …

and

Basilevsky then displays equations for the variables, using Vc not Ac:

What he has written is

i.e.

,

where the Zu are unstandardized, in contrast to Zc and Zs which are standardized. (The Zu were computed from Vc, not from Ac.) I remark that Basilevsky did not compute Ac.

From the correlation matrix:

Now let’s do the correlation matrix, which will give us the third model. Here are the eigenvalues…

and eigenvectors…

Let’s stop there for a moment.

I want to compare the two V’s: recall Vc…

He says, “… even though is positive as a covariance loading, it becomes highly negative as a correlations loadings [sic].”

He correctly points out that the 2nd column of Vc and the second column of Vs would have the same signs if one were multiplied by negative 1, and correctly concludes that we can take care of the sign of by changing the sign of the second column. (At least, I hope that what he said. And is his notation for an element of what i’m calling V.)

Unfortunately, he missed the fact that the 3rd columns cannot be made to have the same signs. The 3rd column of Vc has 2 negative signs; the 3rd column of Vs has only 1.

That is, even though is positive, is negative, when all other signs are the same.

This serves as a reminder that Vc and Vs, eigenvectors of the covariance matrix and of the correlation matrix, do not have a simple relationship to each other.

i will also point out that both Ac and Ar will be derived from Vc by scaling by positive numbers, and As will be derived from Vs by scaling by positive numbers, so this sign discrepancy between the (2,3) elements of Vs and Vc will propagate to the (2,3) elements of Ac and Ar versus As.

Now, lets get back to doing our thing. Get …

and As…

Yes and yes, Basilevsky and I agree. Even to all his signs.

We use As and get the equations

These Xs are standardized, and so are the Zs, in contrast to the Xc and Zu of the first set of equations:

Basilevsky’s model:

Basilevsky has shown us that if we define

,

then we would have the second model,

.

here’s Ar…

And here’s the corresponding equations…

followed immediately by the second set computed earlier…

There are some satisfying similarities. The first columns differ significantly only in the last entry, and even that is .2 vs .36 . For the second column, as we said, if we had changed its sign, all the signs would match up, and in particular the single large terms would then be +0.9324 and +0.9575 .

I should emphasize that he calls both Ar and As, “correlation loadings”, at one point or another.

We could also severely round both As and Ar to confirm their similarity:

As he said, we shoulda changed the sign of the second column of Vc.

Without data, there isn’t a whole lot more we can do here. But I do want to show you one last thing: let’s take a closer look at Ar. We computed it from

and the rows of Ar are supposed to be of unit length. As usual, the fastest way to compute the lengths (more precisely, the squared lengths) of vectors in an array A, is to compute either or . The former gives us dot products of all the rows with each other, including with themselves, while the latter gives us all the dot products among the columns.

So if we compute , the diagonal elements should be 1. We see that they are:

In fact, we see that is equal to our old friend the correlation matrix r:

Don’t be distracted, however, by his scaling the rows. The new basis – which gives us the PCs Zc from Xs – is specified by the columns of Ar, just as the columns of Ac and As were the bases that gave us Zc from Xc, and Zs from Xs. (As I said, all three models are of the same form: . And as I’ve said before, we should keep our eyes on the linear algebra.)

So what does the basis Ar look like? (I’ve said some terrible things about it; let me justify them.) For this, we compute :

That the diagonals are not 1 says the columns are not of unit length; that the off-diagonal elements are not zero says the columns are not mutually orthogonal.

Ar is not a very nice basis at all.

I think it has some redeeming social value, but we’ll look at that when we have data.

Looking at eigenvectors

Let’s check one more thing. The columns of Vc are eigenvectors of the covariance matrix: we should have

.

Do we? We do. (When I compute the difference, I get zero, to within .)

And columns of Ac are eigenvectors of the covariance matrix, too:

.

Similarly, columns of Vs and As are eigenvectors of r, and we should have both

and

.

We do.

(This is why I’m so sure my answers are right: I’ve checked the eigendecompositions.)

What about Ar? Even supposing that its columns are eigenvectors of c, we don’t know what the eigenvalues are. No problem, just see if

r Ar

is proportional to Ar. (That is, see if columns are proportional; see below.)

We have

and when we divide, element by element, by Ar, we get

so we conclude that the columns of Ar are not eigenvectors of r. (Put that the other way: if columns of Ar are eigenvectors, then applying r to them should yield columns which are proportional to the eigenvectors, i.e. proportional to Ar.)

Are the columns of Ar eigenvectors of c? here’s c Ar…

divided by Ar…

Again, no; the columns of Ar are not eigenvectors of c.

Oh, do you need to see how that calculation works when it comes out right? here’s c Ac…

and dividing by Ac gives us…

which says the 4 eigenvalues for the 4 columns are

and that’s exactly what we got a long time ago:

The purpose of this post was to demonstrate Basilevsky’s computations. Next time, I’ll do these kinds of calculations with data.

PCA / FA Basilevsky on standardizing. Discussion.

Introduction and review

Basilevsky presents an extremely interesting idea. For all I know, it’s become common in the last 10-20 years, but I hadn’t seen it in any of the other books we’ve looked at.

I’ll tell you up front that he’s going to normalize the rows of an A matrix, specifically the A matrix computed from the eigendecomposition of the covariance matrix.

I’ll also tell you up front that I don’t see any good reason for doing it, but I’m not averse to finding such a reason someday.

Suppose we have centered data X with variables in columns, and N rows (observations). The covariance matrix is ), and we find its eigendecomposition. We get a matrix V having the eigenvectors as columns, and we construct a diagonal matrix whose entries are the square roots of the eigenvalues.

That is, we have the decomposition

.

Then we define the matrix : each column is an eigenvector, of length .

For the next few posts, I will use a notation of Basilevsky (and Jolliffe): the principal components will be denoted by Z. (Yes, Harman and I almost invariably use , but I need to abandon that notation for a while. Sorry.)

This coincides with Jolliffe’s notation, which was to write

,

with Z as the principal components, X the data, and V the orthogonal eigenvector matrix.

As we have observed many times, the starting point for making sense of that is the change-of-basis equation for a vector:

,

where and are column vectors, x is the old components, z is the new components, and V is the transition matrix. (Its columns are the old components of the new basis vectors. Using the symbol V is not an accident: the matrix V from the eigendecomposition is a transition matrix.)

Applied to the entire matrix X, the change-of-basis equation is

Transpose:

Solve for Z:

But V is orthogonal, , so :

.

Among the very first things we learned in PCA (in Harman) were two: first, that the Z’s we just defined have maximally redistributed the variance of the X’s: the 1st Z variable has the largest variance of any linear combination of the X’s; the 2nd Z variable has the 2nd largest variance of any other linear combination of the X’s, etc. (All subject to the constraint that the total variance of the Z’s is equal to the total variance of the X’s.)

Second, we learned that if we used the transition matrix , then the Z’s defined by it would be both standardized and uncorrelated.

In our present notation (our is Harman’s Z, and our is Harman’s F), we would write Harman’s model as

which, column by column, is precisely the change-of-basis equation

,

with transition matrix A instead of V.

Or we could write

.

(I will try not to mix these up.)

Let me emphasize that the two Z’s in

and

are different. In what follows, I will distinguish them by subscripts (Zu for the first, Zc for the second.)

(It is true that the second set of Z’s could have been found from the first set: just standardize the first set! I had never thought to ask that or check it before, but it would be rather annoying if it were not true. See below.)

We will, in fact, be dealing with both centered data Xc and the corresponding standardized data Xs; we will have corresponding orthogonal eigenvector matrices Vc and Vs; corresponding A matrices Ac and As, and a third one Ar; we will have 3 Z matrices, although a fourth exists in principle.

Zu will be the unstandardized Z from Vc:

Zc will be the standardized Z from Ac:

Zs will be the standardized Z from As: .

(I have tried to minimize the occurence of Z’s as opposed to Zs, and X’s versus Xs, but I can’t reduce the number to zero. Zs and Xs are specific matrices, as opposed to collections of Z’s and X’s, which have apostrophes)

A new basis

We will show that

,

where Ar can be computed from Ac.

That last is a mix of apples and oranges: standardized Xs, but standardized Zc from the covariance matrix, instead of either standardized Xs and standardized Zs from the correlation matrix, or centered Xc and standardized Zc from the covariance matrix.

As I intimated at the beginning, the Ar is Ac with its rows normalized to length 1.

Let’s work that out. The key is that each row of Ac has a sample std deviation as its length.

Let’s see that. To compute the mutual dot products of rows of Ac, we could compute

.

The diagonals of the result would be the squared lengths of rows of Ac (the dot products of rows with themselves).

OK. Let me drop the subscript c for the moment:

,

where c is the covariance matrix.

(the eigendecomposition was .)

So the diagonal elements of are not only the squared lengths of the rows, but also the sample variances of Xc.

We want to scale each row of Ac by its length. To scale columns of V, we post-multiplied by the diagonal matrix . To scale rows of Ac, we will premultiply by a diagonal matrix of inverse standard deviations.

So let = diagonal matrix of square roots of variances, and let be its inverse. (Ignore the possibility of zero variances; after all, we do that every time we talk about the correlation matrix.)

We want to scale the rows of Ac to unit length to get Ar, so we define

.

We have the relationship between Xc and Zc:

.

We standardize Xc by scaling columns by standard deviations, so

.

FInally, transposing

,

QED.

Let me lay out all three equations again (I really like to see them all together):

.

That’s how I would present Basilevsky’s work. He never introduced Ac, and he computed Zc from the unstandardized Zu:

Standardizing Zu

Yes, as I said, the standardized Zc can be computed by standardizing the Zu; and the variances of the Zu are given by the eigenvalues of the covariance matrix . It’s easy enough to confirm.

From our usual definition of Zc,

, so ,

and the old definition of Zu

,

we could write

i.e.

,

QED.

Summary and query

Overall, then, I learned two things. One, that standardized Zc and unstandardized Zu, defined by

and

were related in the sensible way; we could get Zc by standardizing Zu.

Two, and more importantly, I learned that we could write

,

where

.

Fine.

Now, why would we want to?

I don’t know.

Basilevsky got there by looking at

and saying he wanted to standardize both Z and X.

But we already know how to do that:

,

with standardized Xs and standardized Zs.

Why doesn’t he do that?

He does.

And he knows and shows that , and , so by using Ar he has gotten standardized Zc which are different from standardized Zs.

I am perfectly happy with Zs derived from Xs. I am perfectly happy with Zc derived from Xc. I have no idea why we would want to use Ar to relate Zc to Xs. It’s true, but why do we care?

OK, Ac is not an orthonormal basis, but neither is Ar, so that’s not an improvement. OK, the columns of Ac do not have unit length, but they are mutually orthogonal and they are eigenvectors of the covariance matrix. In contrast, the columns of Ar are neither mutually orthogonal , nor are they eigenvectors (of Xc or of Xs).

I find myself wondering if there is some reason for liking this non-orthogonal basis Ar. Does its failure to be orthogonal tell us something? But it didn’t give us anything new. It says (once again) that

,

but Zc is found from Xc and the eigenvectors of the covariance matrix:

,

and Xs is found by standardizing Xc. We don’t need Ar to get the two things it relates.

I think Ar is cute, but I need some reason for going to a basis which isn’t at least orthogonal. Heck, I need a reason for going from the orthogonal Vc to a basis (Ac) whose vectors are not of unit length, but the reason is that it provides standardized Zc.

Finally, I need some reason for abandoning the eigenvectors. I’m in no hurry to substitute Ar for either Ac or As. The eigenvectors have significant properties; what does Ar have to offer that makes up for losing those properties?

If anyone knows….

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.