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 example 4: davis. Summary.

Considering how much I learned from Davis’ little example, I almost hate to summarize it. The following just does not do justice to it. We’ll see if the passage of time – coming back to Davis with fresh eyes after doing the next example – improves the summary.

OTOH, an awful lot of what I got out of Davis came from bringing linear algebra to it, and just plain mucking about with the matrices. And I discovered, having picked up the two chemistry books since I drafted this post, that they made a lot more sense than when I first looked at them. I have learned a lot from working with Davis’ example.

Here’s what I have. My own strong preference, but I’d call it a personal choice, is to use the SVD of the data matrix X:

and compute the A’s and S’s (R-mode and Q-mode loadings and scores):

.

.

.

.

Those equations are worth having as a conceptual basis even if one chooses to do eigendecompositions of or instead of the SVD.

• I would construct a reciprocal basis for .
• And, whenever I construct any basis, I would compute the data wrt that basis.
• And I would probably compute the variances of any new forms of the data.

The Q-mode loadings are the new data wrt the orthogonal eigenvector matrix v (properly named, the columns of v are the right principal vectors of X, but i think of them as the eigenvectors of ). The R-mode loadings (cut down a little) can be thought of as a non-orthonormal basis, and the R-mode scores are the new data wrt the reciprocal basis for .

A lot of what I did with Davis’ example was look at alternative forms and the duality between R-mode and Q-mode. By using the SVD, I choose one particular form, and I have the duality.

Bear in mind that the SVD does not require that X be centered. But that’s a whole ‘nother story.

PCA / FA example 4: davis. Davis & Jolliffe.

What tools did jolliffe give us (that i could confirm)?

1. Z = X A, A orthogonal, X old data variables in columns
2. serious rounding
3. plot data wrt first two PCs
4. how many variables to keep?

but jolliffe would have used the correlation matrix. Before we do anything else, let’s get the correlation matrix for davis’ example. Recall the centered data X:

i compute its correlation matrix:

Now get an eigendecomposition of the correlation matrix; i called the eigenvalues and the eigenvector matrix A. Here’s A:

If we compute and get an identity matrix, then A is orthogonal. (it is.)

the eigenvalues are:

That the third eigenvalue is zero would tell us, if we didn’t already know it from davis, that the data, X, is 2D.

now we can get (1) the new data Z, wrt the orthogonal eigenvector matrix A; it is given by the projection of X onto A (because A is orthogonal, i.e. the basis vectors are orthonormal):

Fascinating.

An understatement. There are 3 nonzero columns.

in sharp contrast to from the centered data X, this new data does not look 2D wrt the orthogonal eigenvector basis. It is still 2D, nothing changed about that, but it isn’t at all clear from the new components of the data. our 3 columns are not linearly independent.

i just gave us all a beautiful lesson: it can be very worthwhile to look at both covariance (or of centered data) and correlation.

BTW, since this is new data, get the variances…

and the total is… 32. we remember that number from davis.

Sure enough, the orthogonal eigenvector matrix A has redistributed the total variance, which is 32 – but it has spread it out over 3 variables instead of 2.

Now we are ready to follow jolliffe. BTW, i will refer to things like “jolliffe’s Z” although jolliffe, of course, never worked out davis’ example! we could try (2) some serious rounding of everything in sight. The correlation matrix: original…

rounded to the nearest .5…

The eigenvector matrix? original…

rounded to the nearest .5…

The new data? original…

Rounded to the nearest .5 …

Rounded to the nearest 1…

i don’t see anything exciting in all that. (but if you’re thinking that i haven’t looked at everything in sight, you’re right. My bad. I’ll make it right. Soon. Very soon.)

As for (3) a plot of the data wrt the first two PCs, that should be: take the first two columns of the new data…

and plot those pairs of points. (this corresponds to davis’ plotting .)

That looks remarkably like a plot of my (mine or davis’ doesn’t matter much, but there are sign differences; nevertheless, the geometry was memorable.) here’s a plot of my :

Clearly we need to look at rounded and Z, and A and v. (i said i’d make it right.) here are jolliffe’s Z …

and my , rounded to the nearest .2 …

The first two columns are not identical, but awfully close, and everything else is identical (to within the usual sign ambiguity).

here are jolliffe’s A …

and my v…

i am surprised at how similar these are. This will not always happen. After all, we saw, in jolliffe when he worked both the correlation and the covariance matrices, that the eigenvectors need not be that similar.

Moving on to (4), how many new variables would jolliffe tell us to keep? recall the eigenvalues of the correlation matrix…

If we keep only those greater than 1 or even those greater than .7, we would keep only one.

Write them as percentages of total…

If we keep those whose cumulative contribution is 70 – 90%, we would keep only one.

a scree graph won’t help at all: we only have two nonzero eigenvalues, so we get a single line segment, and no “elbow” is possible.

What about the broken stick model? recall that if a stick of length 1 is broken into p pieces, then the expected length of the kth largest piece is given by

The expected lengths for a unit stick broken into 3 pieces are…

i.e.

and our fractional eigenvalues were

So, once more, jolliffe would keep only one eigenvector (.8 > .61 but .2 < .278).

Jolliffe notwithstanding, i would keep both nonzero eigenvalues. it just doesn’t seem excessive.

PCA / FA example 4: davis. Reciprocal basis 4.

(this has nothing to do with the covariance matrix of X; we’re back with for the SVD of X.)

in the course of computing the reciprocal basis for my cut down

i came up with the following matrix:

now, is very well behaved. Its two columns are orthogonal to each other:

And each column of is orthogonal to one of the vectors…

and each column of has dot product 1 with the other of the vectors:

i know perfectly well that i could have simply said that is diagonal and that , for my cut-down .

I wanted to write out the orthogonality relations for explicitly. you see, there’s just one tiny problem: i know that is wrong. It is not the reciprocal basis for .

What could be wrong? how can it be so well-behaved but wrong? on a hunch, i computed the normals to the planes spanned by the columns of B and by the columns of , and by the (first two) columns of . in 3D we just use the cross product of two vectors to find their normal.

The cross product of the first two columns of is:

That’s a normal to the plane spanned by the first two columns of . it’s also (anti-)parallel to the third column of v:

.

Of course: the eigenvector basis described by v is orthogonal, so the third one is normal to the plane spanned by the first two.

Similarly, the cross product of the first two columns of B is:

.

(how convenient that the normal vectors have all components equal!) that, too, is (anti-)parallel to the third column of v, so it’s also normal to the plane spanned by the first two columns of .

and the cross product of the columns of ?

.

well, shiver my timbers. that is not like the other two cross products, hence not normal to the plane spanned by the first two columns of .

In contrast to the two vectors of B, the two vectors of do not lie in the same plane as the vectors. They do not span the same plane. They are not a reciprocal basis for the same plane.

is very nice 2D basis. it just happens to span a plane other than the one spanned by . it’s not normal to the third column of v.

It isn’t sufficient to require ; we also need to require that the columns of span the same space – in this case, the same plane – as the columns of !

And how would we realize that our alleged reciprocal basis wasn’t one? well, the projections of X onto the reciprocal basis should be the . if they’re not, something’s wrong.

(guess what? now you know how i discovered that was wrong.)

i would hate weighted eigenvector matrices when there are eigenvalues equal to zero, except that i think i can handle them. this is when i decided to replace zero eigenvalues by 1s, take the inverse transpose, and then take a sub-basis (i.e. the appropriate number of columns).

My own list of things to do now includes constructing the reciprocal basis in precisely that fashion. i choose to construct it because it corresponds to , and let’s me compute the data values for the basis.

PCA / FA example 4: davis. Davis & harman 3.

Hey, since we have the reciprocal basis, we can project onto it to get the components wrt the basis. After all that work to show that the are the components wrt the reciprocal basis, we ought to find the components wrt the basis. And now we know that that’s just a projection onto the reciprocal basis: we want to compute X B.

Recall X:

Recall B:

The product is:

What are the column variances?

Does that surprise you? the 1/3 comes from our using eigenvalues of instead of eigenvalues of the covariance matrix. The new variables have a common variance, but it isn’t 1. this corresponds to our initial work in harman, when we discovered that using the weighted eigenvector matrix gave us new variables with variance 1.

at first glance, it seems not worthwhile to work that out using the covariance matrix, but i suppose it might be worthwhile to come at the covariance matrix from the SVD. (oh, yes! it is!)

Let’s do it.

The covariance matrix of X is

We should just replace X by and find the SVD of X2. but there’s a catch.

Do it anyway. here’s our X2:

Compute the variances…

and their sum … 10.6667 .

Hang on. The data X2 has 1/3 the variance of the original data. We’re trying to find new forms of X (not of X2) with variance 1. but let’s just keep going.

Get the SVD of X2:

here are new u…

new v…

and new w…

Then i compute new A’s and S’s from

.

.

.

.

(using X2, u2, v2, and w2, of course)

are the components of the data X2 wrt the orthogonal eigenvector matrix v2.

Compute the variances…

and their sum is 10.6667 . Yes, the new orthogonal eigenvector matrix v has redistributed the total variance of X2 among two new variables.

What about the R-mode scores?

Compute the variances:

i say again, yikes!

For the record, here are the new

and

From the orthonormal basis v2 , i construct a reciprocal basis for the -weighted matrix . i’ll just use the diagonal of w2 (instead of the since i haven’t computed them) to scale two of the vectors. Recall the new w matrix:

i construct a diagonal matrix using 1/w when possible, 1 otherwise.

then i premultiply that by the new v, and keep only the first two columns; the new reciprocal basis is:

Having the reciprocal basis, we project the data onto it get the new data wrt the basis. But we project the original data. That’s the catch. As i said, it’s the variance of X we’re trying to scale to 1, not the variance of X2. we compute X B2:

And the variances?

.

Indeed. this corresponds even more closely to what we first saw in harman. (if instead we project X2 onto B2, we get a common variance of 1/3 again.)

PCA / FA example 4: davis. Reciprocal basis 2 & 3.

The reciprocal basis for A^R using explicit bases.

Here I go again, realizing that I’m being sloppy. i call

the second data vector, but of course, those are the components of the second data vector. a lot of us blur this distinction a lot of the time, between a vector and its components. So long as we work only with components, this isn’t an issue, but we’re about to write vectors. i wouldn’t go so far as to say that the whole point of matrix algebra is to let us blur that distinction, but it’s certainly a major reason why we do matrix algebra. but for now, let’s look at the linear algebra, distinquishing between vectors and their components.

i need a name for the second data vector; let’s call it s. we will write it two ways, with respect to two bases, and show that the two ways are equivalent.

We have the original basis. We never write it, because the basis vectors have components (1,0,0), (0,1,0), and (0,0,1). let me call this basis . because it’s an orthonormal basis, it is its own reciprocal basis, and i will call the reciprocal basis . yes, i use the same symbol e, but the original basis has subscripts, the reciprocal basis has superscripts.

We have the basis described by . let me call those vectors . we found the reciprocal basis B, and i will call it . same convention: subscripts or superscripts on the common symbol f.

let us recall the (transition matrix B) for the reciprocal basis:

let us also recall the design matrix X:

When we say the components of the second data vector are

we are asserting that the second data vector is

.

When we say that the first reciprocal basis vector has components given by the first column of B…

we are asserting that the first reciprocal basis vector is

Similarly, the second reciprocal basis vector, whose components are the second column of B, is

Let us recall (the first two columns of) :

When we say that the second row of …

are the components of the second data vector wrt the reciprocal basis , we are saying that

.

but we defined the vector s as

.

The two forms of s should be the same vector. We start with…

.

Plug in the expressions for and …

and Mathematica® simplifies that to

,

just what it should be.

The reciprocal basis for A^R using tensors.

Third, i can write this another way, and some of you will have seen it this way. If you have, i just want to refresh your memory. Using tensor notation (and the einstein summation convention), we could have written

and .

which is another way of saying that the components of s wrt the reciprocal basis are and the components of s wrt the basis are . It is absolutely crucial that corresponding components and vectors do not both have superscripts or subscripts, but are mixed.

what is significant are these equations

and

,

which I offer as reminders; i’m not going to work them out for you at this time. they are another way of writing what we’ve been saying all along: dotting a vector into some basis vectors gives you the components wrt the reciprocal basis.

PCA / FA example 4: davis. Scores & reciprocal basis 1.

The reciprocal basis for , using matrix multiplication.

i keep emphasizing that is the new data wrt the orthogonal eigenvector matrix; and i hope i’ve said that the R-mode scores are not the new data wrt the weighted eigenvector matrix . In a very real sense, the R-mode scores do not correspond to the R-mode loadings . This is important enough that i want to work out the linear algebra. In fact, i’ll work it out thrice. (And i’m going to do it differently from the original demonstration that is the new data wrt v.)

Of course, people will expect us compute the scores , and I’ll be perfectly happy to. I just won’t let anyone tell me they correspond to the .

First off, we have been viewing the original data as vectors wrt an orthonormal basis. Recall the design matrix:

Let’s take just the second observation:

Those are the components of a vector.

Second, we have an orthonormal basis defined by the orthogonal eigenvector matrix v. Recall :

for the second row, in particular, we have new components….

We’ve already seen that old and new components (x and y, resp) are related by

,

although we wrote the whole thing as . yes? we confirm it simply by computing…

x

=

which is .

The second row of represents the same vector as the second row of X, but the components are different because the basis vectors are different.

All that was review, albeit from a slightly different point of view. let’s move on to the weighted eigenvector matrix ; it is not an orthonormal basis because the basis vectors are not of length 1, although they are orthogonal to each other. Here’s my from the SVD:

i just want the first two columns. Here’s my cut down…

and my cut down:

I’ve been implying, if not saying, that the projections of X onto , namely (which are the R-mode scores ), are the components of the data wrt a basis other than . it’s what may best be called the reciprocal basis.

(if you haven’t looked at “the reciprocal basis” in http://rip94550.wordpress.com/2008/03/21/transpose-matrix-adjoint-operator-2/, you might want to.)

Briefly, a reciprocal basis is constructed to make the dot product of vectors work out right when we have a non-orthonormal basis. What we do is get the components of one vector wrt the basis, the other vector wrt the reciprocal basis, and then we can compute their dot product in the usual way, as the sum of products of corresponding entries; i.e. as a row vector times a column vector.

If we have a transition matrix A, then the transition matrix B for the reciprocal basis is B = . (because it gives us , which says that the dot product of a row of with a column of A is either 0 or 1. it’s equivalent to saying that the dot product of a column of B with a column of A is either 0 or 1. either way, the column vectors of B and of A are pairwise orthogonal.)

we have a problem, however: we have a matrix which is not invertible. Another damned dragon to cope with.

We could look for a matrix B such that . unfortunately, such matrices are not unique.

Let’s construct the reciprocal basis geometrically. As was constructed from v by multiplying columns by we will construct B by dividing v by , but we’ll just do it for the first two columns. The quickest way to accomplish that is to post-multiply by a diagonal matrix whose entries are , when possible and 1 otherwise… the diagonal eigenvalue matrix is…

so i form a new diagonal matrix

Note the 1 instead of 0 in the (3, 3) slot.

recall v…

and compute the product of v and the new diagonal matrix

x

=

and since i only want the first two columns, i cut that down to

and call it B.

It is true that . (so the columns of B are mutually orthogonal to the columns of A.)

Now, B is the transition matrix for the reciprocal basis of . was defined as , the projections of the data onto the basis. i claim that is the components of the data wrt B, the reciprocal basis for .

i want to confirm this just for the second data vector, namely that

,

i.e. the old components are found by applying the transition matrix B to the new components. the RHS is

x

=

and, once again, we have found the old components of the second observation. More importantly, we have exhibited a basis wrt which the R-mode scores represent the new components of the data.

PCA / FA example 4: davis. Davis & harman 2.

What about redistributing the variance?

We learned from harman that if we use the orthogonal eigenvector matrix v to compute new data, the result has redistributed the original variance according to the eigenvalues; and we also learned that if we use the weighted eigenvector matrix instead, we get new variables of unit variance.

(i think i thought this required us to start with the correlation matrix; it does not.)

let’s see this again for davis’ centered data X:

We have variances…

or…

and the sum of the variances is…

32.

Now is a good time to recall :

The sum of the diagonal elements is 96, and if we divide by 3 (= N-1, for the covariance matrix instead of for ), gee, we get 32.

stop for a moment. We have two, admittedly proportional, sets of eigenvalues, one set from the covariance matrix, the other from . this gives us two different possible sets of weights to apply to the eigenvector matrix v. (recall that v is an eigenvector matrix both for and for the covariance matrix of X, as well as coming out of the SVD of X.)

Let’s stay with what i know: our new data, using the orthogonal eigenvector matrix v, is .

The variances are…

and we don’t even need to consciously add those: the sum is 32. We have redistributed the original variance among our two new variables.

Now, everyone seems to say that the -weighted eigenvectors have variances equal to the eigenvalues. What those eigenvectors really have is lengths equal to their .

here indeed be dragons.

as linear algebra is my touchstone in general for all of this, that redistributed variance using the orthogonal eigenvector matrix is my vorpal blade (dragon, jabberwock, no big deal – snicker-snack either way). any different choice of basis must lead to different results.

we looked at that redistribution of variance in harman. We learned then that the correct way to get new data with redistributed variance is to use the orthogonal eigenvector matrix. We just reconfirmed it for davis’ example.

Instead of actually kicking one of the dragons awake, let’s try counting its scales while it sleeps. If the R-mode scores are some form of data (and they are), what are their variances?

my is

the variances of the first two columns are

yikes! if nothing else, that should convince you that the R-mode scores are weird.Life is so much simpler if we use the orthogonal eigenvector matrix. The new data is trivial to compute (X v), and we can see that it has redistributed the variance of the original data.Trying to use any weighted eigenvector matrix is dicey if there is a zero eigenvalue (because our transition matrix is not invertible); and confusing if we used instead of the covariance matrix (because they have different eigenvalues).i will eventually calculate the data wrt both of the weighted eigenvector matrices. The dragons will be awake and aloft and pretty to look at.

PCA / FA example 4: davis. Davis & harman 1.

Let us now recall what we saw harman do. Does he have anything to add to what we saw davis do? if we’re building a toolbox for PCA / FA, we need to see if harman gave us any additional tools.

here’s what harman did or got us to do:

• a plot of old variables in terms of the first two new ones
• redistribute the variance of the original variables
• his model, written Z = A F, was (old data, variables in rows) = (some eigenvector matrix) x (new data)

We should remind ourselves, however, that harman used an eigendecomposition of the correlation matrix; and the eigenvectors of the correlation matrix are not linearly related to the eigenvectors of the covariance matrix (or, equivalently, of .) having said that, i’m going to go ahead and apply harman’s tools to the davis results. i will work from the correlation matrix when we get to jolliffe.

What exactly was that graph of harman’s? he took the first two columns of the weighted eigenvector matrix, then plotted each row as a point. He had 5 original variables. His plot

shows that the old 1st and 3rd variables, for example, are defined very similarly in terms of the first two new variables. let me emphasize that this is a description of defining variables, not a display of the resulting data. Simple counting can help with the distinction, e.g. 3 variables but 4 observations.

from davis we only have 3 variables to start, and our weighted eigenvector matrix has only two nonzero columns; and it’s one form of what we called . we have to be careful, because there are sign differences among the various forms of ; not because the definitions are different, but because eigenvectors are not unique. In particular, my u and v from the SVD have some different signs from davis’ U and V from eigendecompositions.

most of the time, i will use my SVD results, but right now let’s compare davis with harman. Davis had as…

Davis showed us a plot of these numbers viewed as three points.

So the R-mode loadings tell us how the 3 old variables are described by the new ones; we see that a plot of corresponds to harman’s plot.

i want to emphasize that even the orthogonal eigenvector matrix shows us that the data is really only 2D. Don’t think that we had to scale the eigenvectors in order to get 2D.

my orthogonal eigenvector matrix v is

but wrt that basis, the new data is :

i.e. the new data is described by the first two eigenvectors, hence it lies in a plane. but then so does the original data, because the transformation between them is linear (a plane cannot be mapped to 3D, although it could be mapped to a line). we could plot the data in 3D and rotate the point of view until we see the plane edge on.

Do davis and harman’s models resemble each other?

Let’s recall the SVD form of the A’s and S’s.

.

.

.

.

Then

.

in particular,

and then .

is harman’s model (Z = A F, where , A is a weighted eigenvector matrix, F is the new data), while is jolliffe’s (Z = X A, where Z is the new data, X is the original data, and A is the orthogonal eigenvector matrix). i think it would be unfair to say that davis is using either of those models, especially since harman and jolliffe (usually) work from the correlation matrix.

i would summarize davis with all four A’s and S’s.

i can’t, however, say too many times that and can be constructed using u or v:

,

and

.

i see no grounds for associating with just u or just v. and no grounds for associating either u or v primarily with R-mode or Q-mode; they’re equally important to the SVD.

i would probably focus on v and : one defines new variables, the other gives me the associated new data.

PCA / FA example 4: Davis review (3)

Let’s talk about centering the data. Recall the questions:

• is there any chance that we can always center both the columns and the rows?
• is there any chance that we can always standardize both the columns and the rows?
• if we can’t have our cake and eat it too, should we give up the duality between R-mode and Q-mode?
• what if rows of the data matrix add up to 1 or 100 (i.e. the variables are percentages)?

(My third question is poorly phrased. we always have the duality between A’s and S’s; what i feared losing was the idea that Q-mode is just R-mode applied to the transpose ; but what if X doesn’t have row-centered data?)

The third question is also a red herring. We can always make both the columns and the rows add up to zero. Proving it isn’t too hard; some convenient notation might simplify things. If the data is

then

and

are good symbols for the column means and row means, respectively. The grand mean (the mean of all the matrix values = mean of column means = mean of row means) would be denoted

BTW, it’s very handy that the mean of column means = mean of row means; we’ll see how handy, soon.

The only subtlety is that what we have here are row and column means of the original data, instead of, for example, column means of the original data and row means of the column-centered data. We can get “doubly-centered data” by computing

Why are we adding the grand mean? Because we subtracted it twice. If we had taken row means of the column centered data, we wouldn’t add the grand mean back in; but i chose, and i’m sure most people choose, to take row means of the original data.

I leave the proof that we can always get doubly-centered data for the interested reader. Really, it shouldn’t be difficult.

It means that we can always do both the R-mode and Q-mode analyses on a doubly-centered design matrix X.

That had bothered while i was looking at Q-mode. I knew that davis’ data was row-centered as well as column-centered (had zero-mean rows as well as zero-mean columns, resp.), but i just assumed he had picked the numbers that way in order to avoid discussing the case when the rows were not zero-mean.

We can always have both row-centered and column centered simultaneously.
(We can have our cake and eat it too, so we don’t have to worry about Q-mode being done with non-centered data.)

NOTE that we cannot do that with standardized data (or, therefore, with the correlation matrices). That is, we cannot standardize both the columns and the rows. The unit column-variances get messed up when we standardize the rows after standardizing the columns.

As I draft this, it is an open question for me whether we can always arrange to have column-centered data with row sums equal to 1. Let me follow some great advice, scary as it is: let me guess. (from John A. Wheeler, the great physics teacher: if my guess turns out right, i have reinforced my intuition; if my guess turns out wrong, i correct my intuition.)

I think we can do it. Adding that grand mean back in looks promising for adding 1 back in.

(I won’t string you along any further: my intuition was wrong. We’ll see what happens instead. It’s almost as good.)

Let’s try an example. I start with a real simple set of numbers:

Square the first row, square-root the second row, leave the third alone:

While i’ve got it in this form, get what will be the row means.

(Much of Mathematica® works on rows; I’ve just been working with .) Transpose. This is my example data matrix.

Get the column means.

The grand mean is any of the following: the mean of row means… the mean of column means… or the mean of all 12 values. No matter how we compute it, we get

6.84664 .

We get doubly-centered data by

which gives us

(Yes, the column means are zero and the row means are zero.)

That was easy enough.

Let’s back up and look at just column-centered data from the original.

The column means, of course, are zero… the row means are not (these, of course, are the row means of the column centered data, not of the original data)…

but the grand mean is zero. (go ahead, compute it if you must.)

The grand mean is zero because it is the mean of the column means, each of which is zero. This implies that the mean of the row means is zero. In particular, since 1 ≠ 0, we cannot have column-centered data with constant row sums = 1. (I said it was handy that the grand mean could be computed more than one way.)

maybe I should emphasize what we just saw: the act of column-centering the data implies that the sum of the new row means is zero. If each row mean is the same, then that common value must be zero. That example did not have constant row sums, but we see what must happen if it had.

I probably do not need to work an example, but i’m going to. And i’m going to do it because I had actually worked this example first; then i realized what had happened, and i told you you up front. So suppose we had started with a matrix where the rows each add up to 1. Go back to the raw data again….

but change the rows. Divide each row by its mean, and by 3. (I don’t want row mean = 1, I want row sum = 1.)

Get the column means…

and having said that we want row sums = 1, we confirm that each row mean is 1/3, because it’s an easier command to Mathematica®:

I then construct column-centered data:

We could confirm that the column means are zero; and we compute the row means…

We automatically got row-centered data, as expected.

We could look at standardized data, but i’ll leave this to you, too: a counterexample should be convincing enough (and I’ve given you an matrix to play with). We can’t preserve unit variance of the columns when we then standardize the rows.

Summary:

1. we can always get doubly-centered data.
2. my guess was wrong: we cannot have column-centered with constant row sums, unless the constant is 0.
3. but if the rows do have a common row mean (equivalently, a common row sum), then just centering the columns automatically gives us doubly-centered data.
4. we cannot always get doubly-standardized data.
5. if i center the columns, then i would be inclined to center the rows: use doubly-centered data in preference to only column-centered (assuming, as i prefer, that the variables are columns.)

Let me elaborate on (5). For a matrix X with variables in columns, column-centered data means that is proportional to the covariance matrix of X. If we want to be proportional to the covariance matrix of , then we need row-centered data. If we want both and to be proportional to covariance matrices, then we need doubly-centered data.