0
$\begingroup$

I'm trying to understand, in simple terms, the following example copied from prcomp in R:

C <- chol(S <- toeplitz(.9 ^ (0:31))) # Cov.matrix and its root
all.equal(S, crossprod(C))
set.seed(17)
X <- matrix(rnorm(32000), 1000, 32)
Z <- X %*% C  ## ==>  cov(Z) ~=  C'C = S
all.equal(cov(Z), S, tol = 0.08)
pZ <- prcomp(Z, tol = 0.1)

I gather from the included comment that the toeplitz matrix S is suitable as covariance matrix, and C is its root. The second line proves this. What I don't understand is why multiplying a matrix of random numbers X by C gives a "data set" that is (somehow more) desirable as a demonstration of prcomp. Why not simply use X to demonstrate prcomp? Logically, it seems to me that Z must have some special "structure" that guarantees something about the PCA results. However, my matrix algebra is limited, and I don't understand what's going on here. Can someone explain in conceptual terms why in the example it is desirable to create Z?

$\endgroup$
6
  • $\begingroup$ The answer is right there in the comment to the code, which asserts the covariance of $C$ is $S.$ By contrast, $X$ is constructed from a distribution with a rather "uninteresting" covariance matrix--a unit diagonal. $\endgroup$
    – whuber
    Commented Mar 16, 2019 at 16:51
  • $\begingroup$ Sorry, I'm confused. Did you mean cov(Z) = S? I can see that is true (within limits), but it doesn't help me understand. On your second sentence: diag(cov(X)) is roughly 1.0, I can see that, and plot(prcomp(X)) shows a slow (uninteresting) decline. But, plot(prcomp(Z)) falls off more like a real data set might, as does plot(prcomp(C)) as well as plot(prcomp(S)). So why not demo with S or C instead of creating Z? $\endgroup$ Commented Mar 16, 2019 at 18:04
  • $\begingroup$ So, looking at how I phrased the question originally, I see that matrices Z, C, S all have some "special structure" by which I mean you get an interesting scree plot. X has that structure too by virtue of being constructed using C. So why all the extra work when plot(prcomp(S)) would have done the trick? I realize I am asking you to explain someone else's thinking here, just trying to get a better understanding. $\endgroup$ Commented Mar 16, 2019 at 18:11
  • $\begingroup$ X has nothing to do with C: it is generated randomly according to a collection of independent standard Normal variables. The purpose of the example is to exhibit the results of a PCA on a plausible collection of data Z that have sizable correlations. prcomp(S) wouldn't do that. $\endgroup$
    – whuber
    Commented Mar 16, 2019 at 18:31
  • 1
    $\begingroup$ OK, I think I get it. While plot(prcomp(S)) and plot(prcomp(Z)) give qualitatively similar scree plots, the magnitude (y-axis) for the scree plot of Z is much higher, and of course we want to explain much more variance. So creating and using Z gives a better example. I may summarize in a self-answer later. Or you can add an answer. Thank you so much! $\endgroup$ Commented Mar 16, 2019 at 19:12

1 Answer 1

1
$\begingroup$

...Answering my own question for the sake of having an answer, with all thanks due to @whuber.

The key question is "what is it about matrix Z that makes it a good demo data set/why go to all the trouble they did to make Z?"

Short answer: it has "structure" and it is of appropriate scale.

Here's a tour to give a fuller answer:

First, X as defined in the example is of course a draw from a normal distribution and does not have an interesting structure. By "interesting structure" I mean it would give a scree plot resembling something you'd like to see in a real data set. Here is the scree plot created via plot(prcomp(X)):

enter image description here

Notice that the explained variance falls off very slowly, and there is no elbow evident. Further, the amount of variance explained is very low. This is on the y-axis. Unfortunately the built-in scree plot gives absolute values for the variance explained. To see the percent variance explained, do summary(prcomp(X)). In this case the PCs are only explaining about 4% each. So X does not give a good or realistic demonstration of PCA.

In the example it appears that the desired structure comes from the construction of the toeplitz matrix. It's easy to see what this is via the example in ?toeplitz (look at the diagonals):

> x <- 1:5
> toeplitz (x)
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    2    3    4    5
[2,]    2    1    2    3    4
[3,]    3    2    1    2    3
[4,]    4    3    2    1    2
[5,]    5    4    3    2    1

One can see there is some structure to this matrix. Do let's do the same steps as we did for X using S, the toeplitz matrix created in ?prcomp. Here's the scree plot which looks like it came from a real data set, in that it drops off quickly.

enter image description here

And inspection of summary(prcomp(S)) shows that three PCs explains > 95% of the variation. So this is not a bad demonstration. Perhaps the author of ?prcomp could have stopped here.

However, they went ahead and found the root of S, C, via Cholesky decomposition, and then post-multiplied X with it to give Z (Z = X %*% C). In the comments they point out that cov(Z) ~= S. That's true (they prove it), I guess it's important because it shows the structure of S is carried through. For instance, plot(prcomp(X %*% S)) gives something visually similar to plot(prcomp(S)) though the vertical scale is quite different (the percent variance explained is roughly similar though).

The only reason I can think of that the authors go all the way to Z in their demo is "noise". The scree plots and summaries of prcomp(S or C or Z) are similar. The difference however is in the range of the scores (note that range on class(prcomp) is the range of the scores, which roughly parallel the range of the raw data):

> range(prcomp(S))
[1] -1.577132  1.577132
> range(prcomp(C))
[1] -1.7389697  0.7116142
> range(prcomp(Z))
[1] -11.28021  13.54936

So I guess that a data set like Z is less noisy and more desirable as a demonstration if all other factors are equal (specifically, the trend in variance explained is similar). If this is applicable to real data sets, I don't know what it is called (and in the realm of real data sets we usually have the data we have, and that's what must be analyzed, there is no other choice).

$\endgroup$
3
  • $\begingroup$ @whuber I'd love to have your feedback on this answer. $\endgroup$ Commented Mar 17, 2019 at 13:56
  • 1
    $\begingroup$ +1 It's a great account of your thought process. The one missing thing seems to be this: in applications, all we have available is the data matrix $Z.$ We don't know $S.$ The code example provides you an advantage you don't ordinarily have: it enables you to compare the PCA of the data to that of the true underlying covariance. This is in keeping with all statistical theory: to figure out how well a procedure works, study how well it works when the underlying distribution is known. $\endgroup$
    – whuber
    Commented Mar 17, 2019 at 16:30
  • $\begingroup$ @whuber Excellent point, which I am now in a position to appreciate. I'll accept my own answer so the case is closed, but with great thanks to you. This has pushed my understanding along nicely. $\endgroup$ Commented Mar 17, 2019 at 16:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.