Take the 2-minute tour ×
Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. It's 100% free, no registration required.

I'm working in R through an excellent PCA tutorial by Lindsay I Smith and am getting stuck in the last stage. The R script below takes us up to the stage (on p.19) where the original data is being reconstructed from the (singular in this case) Principal Component, which should yield a straight line plot along the PCA1 axis (given that the data only has 2 dimensions, the second of which is being intentionally dropped).

d = data.frame(x=c(2.5,0.5,2.2,1.9,3.1,2.3,2.0,1.0,1.5,1.1),
               y=c(2.4,0.7,2.9,2.2,3.0,2.7,1.6,1.1,1.6,0.9))

# mean-adjusted values 
d$x_adj = d$x - mean(d$x)
    d$y_adj = d$y - mean(d$y)

# calculate covariance matrix and eigenvectors/values
(cm = cov(d[,1:2]))

#### outputs #############
#          x         y
# x 0.6165556 0.6154444
# y 0.6154444 0.7165556
##########################

(e = eigen(cm))

##### outputs ##############
# $values
    # [1] 1.2840277 0.0490834
    #
    # $vectors
#          [,1]       [,2]
# [1,] 0.6778734 -0.7351787
# [2,] 0.7351787  0.6778734
###########################


# principal component vector slopes
s1 = e$vectors[1,1] / e$vectors[2,1] # PC1
s2 = e$vectors[1,2] / e$vectors[2,2] # PC2

plot(d$x_adj, d$y_adj, asp=T, pch=16, xlab='x', ylab='y')
abline(a=0, b=s1, col='red')
abline(a=0, b=s2)

enter image description here

# PCA data = rowFeatureVector (transposed eigenvectors) * RowDataAdjust (mean adjusted, also transposed)
feat_vec = t(e$vectors)
row_data_adj = t(d[,3:4])
final_data = data.frame(t(feat_vec %*% row_data_adj)) # ?matmult for details
names(final_data) = c('x','y')

#### outputs ###############
# final_data
#              x           y
# 1   0.82797019 -0.17511531
# 2  -1.77758033  0.14285723
# 3   0.99219749  0.38437499
# 4   0.27421042  0.13041721
# 5   1.67580142 -0.20949846
# 6   0.91294910  0.17528244
# 7  -0.09910944 -0.34982470
# 8  -1.14457216  0.04641726
# 9  -0.43804614  0.01776463
# 10 -1.22382056 -0.16267529
############################

# final_data[[1]] = -final_data[[1]] # for some reason the x-axis data is negative the tutorial's result

plot(final_data, asp=T, xlab='PCA 1', ylab='PCA 2', pch=16)

enter image description here

This is as far as I've got, and all OK so far. But I can't figure out how the data is obtained for the final plot - the variance attributable to PCA 1 - which Smith plots as:

enter image description here

This is what I've tried (which ignores adding the original means):

trans_data = final_data
trans_data[,2] = 0
row_orig_data = t(t(feat_vec[1,]) %*% t(trans_data))
plot(row_orig_data, asp=T, pch=16)

.. and got an erronous:

enter image description here

.. because I've lost a data dimension somehow in the matrix multiplication. I'd be very grateful for an idea what's going wrong here.


* Edit *

I wonder if this is the right formula:

row_orig_data = t(t(feat_vec) %*% t(trans_data))
plot(row_orig_data, asp=T, pch=16, cex=.5)
abline(a=0, b=s1, col='red')

But I'm a little confused if so because (a) I understand the rowVectorFeature needs to be reduced to the desired dimensionality (the eigenvector for PCA1), and (b) it doesn't line up with the PCA1 abline:

enter image description here

Any views much appreciated.

share|improve this question

3 Answers 3

up vote 4 down vote accepted

You were very very nearly there and got caught by a subtle issue in working with matrices in R. I worked through from your final_data and got the correct results independently. Then I had a closer look at your code. To cut a long story short, where you wrote

row_orig_data = t(t(feat_vec[1,]) %*% t(trans_data))

you would have been ok if you had written

row_orig_data = t(t(feat_vec) %*% t(trans_data))

instead (because you'd zeroed out the part of trans_data that was projected on the second eigenvector). As it was you were trying to multiply a $2\times1$ matrix by a $2\times10$ matrix but R did not give you an error. The problem is that t(feat_vec[1,]) is treated as $1\times2$. Trying row_orig_data = t(as.matrix(feat_vec[1,],ncol=1,nrow=2) %*% t(trans_data)) would have given you a non-conformable arguments error. The following, possibly more along the lines of what you intended, would also have worked

row_orig_data = t(as.matrix(feat_vec[1,],ncol=1,nrow=2) %*% t(trans_data)[1,])

as it multiplies a $2\times1$ matrix by a $1\times10$ matrix (note that you could have used the original final_data matrix here). It isn't necessary to do it this way, but it's nicer mathematically because it shows that you are getting $20=2\times10$ values in row_orig_data from $12=2\times1 + 1\times10$ values on the right hand side.

I've left my original answer below, as someone might find it useful, and it does demonstrate getting the required plots. It also shows that the code can be a bit simpler by getting rid of some unnecessary transposes: $(XY)^T=Y^TX^T$ so t(t(p) %*% t(q)) = q %*% t.

Re your edit, I've added the principal component line in green to my plot below. In your question you got had the slope as $x/y$ not $y/x$.


Write

d_in_new_basis = as.matrix(final_data)

then to get your data back in its original basis you need

d_in_original_basis = d_in_new_basis %*% feat_vec

You can zero out the parts of your data that are projected along the second component using

d_in_new_basis_approx = d_in_new_basis
d_in_new_basis_approx[,2] = 0

and you can then transform as before

d_in_original_basis_approx = d_in_new_basis_approx %*% feat_vec

Plotting these on the same plot, together with the principal component line in green, shows you how the approximation worked.

plot(x=d_in_original_basis[,1]+mean(d$x),
         y=d_in_original_basis[,2]+mean(d$y),
     pch=16, xlab="x", ylab="y", xlim=c(0,3.5),ylim=c(0,3.5),
     main="black=original data\nred=original data restored using only a single eigenvector")
points(x=d_in_original_basis_approx[,1]+mean(d$x),
           y=d_in_original_basis_approx[,2]+mean(d$y),
       pch=16,col="red")
points(x=c(mean(d$x)-e$vectors[1,1]*10,mean(d$x)+e$vectors[1,1]*10), c(y=mean(d$y)-e$vectors[2,1]*10,mean(d$y)+e$vectors[2,1]*10), type="l",col="green")

enter image description here

Let's rewind to what you had. This line was ok

final_data = data.frame(t(feat_vec %*% row_data_adj))

The crucial bit here is feat_vec %*% row_data_adj which is equivalent to $Y=S^TX$ where $S$ is the matrix of eigenvectors and $X$ is your data matrix with your data in rows, and $Y$ is the data in the new basis. What this is saying is that the first row of $Y$ is the sum of (rows of $X$ weighted by the first eigenvector). And the second row of $Y$ is the sum of (rows of $X$ weighted by the second eigenvector).

Then you had

trans_data = final_data
trans_data[,2] = 0

This is ok: you're just zeroing out the parts of your data that are projected along the second component. Where it goes wrong is

row_orig_data = t(t(feat_vec[1,]) %*% t(trans_data))

Writing $\hat Y$ for the matrix of data $Y$ in the new basis, with zeros in the second row, and writing $\mathbf{e}_1$ for the first eigenvector, the business end of this code t(feat_vec[1,]) %*% t(trans_data) comes down to $\mathbf{e}_1 \hat Y$.

As explained above (this is where I realised the subtle R problem and wrote the first part of my answer), mathematically you are trying to multiply a $2\times1$ vector by a $2\times10$ matrix. This doesn't work mathematically. What you should do is take the first row of $\hat Y$ = the first row of $Y$: call this $\mathbf{y}_1$. Then multiply $\mathbf{e}_1$ and $\mathbf{y}_1$ together. The $i$th column of the result $\mathbf{e}_1\mathbf{y}_1$ is the eigenvector $\mathbf{e}_1$ weighted by the 1st coordinate only of the $i$th point in the new basis, which is what you want.

share|improve this answer
    
Thanks TooTone this is very comprehensive, and resolves the ambiguities in my understanding of the matrix calculation and role of featureVector in the final stage. –  geotheory Mar 18 at 15:36
    
Great :). I answered this question because I am studying the theory of SVD / PCA at the moment and wanted to get to grips with how it works with an example: your question was good timing. After working through all the matrix calculations I was a bit surprised that it turned out to be an R problem -- so I'm glad you appreciated the matrices aspect of it too. –  TooTone Mar 18 at 15:39

I think you have the right idea but stumbled over a nasty feature of R. Here again the relevant code piece as you've stated it:

trans_data = final_data
trans_data[,2] = 0
row_orig_data = t(t(feat_vec[1,]) %*% t(trans_data))
plot(row_orig_data, asp=T, pch=16)

Essentially final_data contains the coordinates of the original points with respect to the coordinate system defined by the eigenvectors of the covariance matrix. To reconstruct the original points one therefore has to multiply each eigenvector with the associated transformed coordinate, e.g.

(1) final_data[1,1]*t(feat_vec[1,] + final_data[1,2]*t(feat_vec[2,])

which would yield the original coordinates of the first point. In your question you set the second component correctly to zero, trans_data[,2] = 0. If you then (as you already edited) calculate

(2) row_orig_data = t(t(feat_vec) %*% t(trans_data))

you calculate formula (1) for all points simultaneously. Your first approach

row_orig_data = t(t(feat_vec[1,]) %*% t(trans_data))

calculates something different and only works because R automatically drops the dimension attribute for feat_vec[1,], so it is not a row vector anymore but treated as a column vector. The subsequent transpose makes it a row vector again and that's the reason why at least the calculation does not produce an error, but if you go through the math you will see that it is something different than (1). In general it is a good idea in matrix multiplications to suppress dropping of the dimension attribute which can be achieved by the drop parameter, e.g. feat_vec[1,,drop=FALSE].

Your edited solution seems correct, but you calculated the slope if PCA1 wrongly. The slope is given by $\Delta y / \Delta x$, hence

s1 = e$vectors[2,1] / e$vectors[1,1] # PC1
s2 = e$vectors[2,2] / e$vectors[1,2] # PC2
share|improve this answer
    
Thanks very much Georg. You're right about the PCA1 slope. Very useful tip also about the drop=F argument. –  geotheory Mar 18 at 15:32

After exploring this exercise you can try the easier ways in R. There are two popular functions for doing PCA: princomp and prcomp. The princomp function does the eigenvalue decomposition as done in your exercise. The prcomp function uses singular value decomposition. Both methods will give the same results almost all of the time: this answer explains the differences in R, whereas this answer explains the math. (Thanks to TooTone for comments now integrated into this post.)

Here we use both to reproduce the exercise in R. First using princomp:

d = data.frame(x=c(2.5,0.5,2.2,1.9,3.1,2.3,2.0,1.0,1.5,1.1), 
               y=c(2.4,0.7,2.9,2.2,3.0,2.7,1.6,1.1,1.6,0.9))

# compute PCs
p = princomp(d,center=TRUE,retx=TRUE)

# use loadings and scores to reproduce with only first PC
loadings = t(p$loadings[,1]) 
scores = p$scores[,1] 

reproduce = scores %*% loadings  + colMeans(d)

# plots
plot(reproduce,pch=3,ylim=c(-1,4),xlim=c(-1,4))
abline(h=0,v=0,lty=3)
mtext("Original data restored using only a single eigenvector",side=3,cex=0.7)

biplot(p)

enter image description here enter image description here

Second using prcomp:

d = data.frame(x=c(2.5,0.5,2.2,1.9,3.1,2.3,2.0,1.0,1.5,1.1), 
               y=c(2.4,0.7,2.9,2.2,3.0,2.7,1.6,1.1,1.6,0.9))

# compute PCs
p = prcomp(d,center=TRUE,retx=TRUE)

# use loadings and scores to reproduce with only first PC
loadings = t(p$rotation[,1])
scores = p$x[,1]

reproduce = scores %*% loadings  + colMeans(d)

# plots
plot(reproduce,pch=3,ylim=c(-1,4),xlim=c(-1,4))
abline(h=0,v=0,lty=3)
mtext("Original data restored using only a single eigenvector",side=3,cex=0.7)

biplot(p)

enter image description here enter image description here

Clearly the signs are flipped but the explanation of variation is equivalent.

share|improve this answer
    
Thanks mrbcuda. Your biplot looks identical to Lindsay Smith's so I presume he/she used the same method 12 years ago! I'm also aware of some other higher level methods, but as you rightly point out this is an exercise to make the underlying PCA maths explicit. –  geotheory Mar 18 at 15:44

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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