I have a bunch of points in $\mathbb{R}^3$ that I would like to translate and rotate so that their center is at the origin and the variance along the $x$ and $y$ axes are maximal (greedy, and in that order). To accomplish this I am trying to use python's principal components analysis algorithm. It is not behaving as I expect it to, most likely due to some misunderstanding about what PCA actually does on my part.

The Problem: When I center and then rotate the data, the variance along the third component is greater than along the second. This means that, once centered and rotated, there is more variance in the data along the $z$ axis than there is along the $y$. In other words, the rotation is not the correct one.

What I am Doing: Python's PCA routine returns an object (say myPCA) with several attributes. myPCA.Y is the data array, but centered, scaled, and rotated (in that order). I do not want the data to be scaled. I simply want a translation and a rotation.

import numpy as np                                                              
from matplotlib.mlab import PCA                                                 

# manufactured data producing the problem                                       
data_raw  = np.array([                                                          
                     [80.0, 50.0, 30.0],                                        
                     [50.0, 90.0, 60.0],                                        
                     [70.0, 20.0, 40.0],                                        
                     [60.0, 30.0, 45.0],                                        
                     [45.0, 60.0, 20.0]                                         

# obtain the PCA                                                                
myPCA = PCA(data_raw)                                                           

# center the raw data                                                           
centered = np.array([point - myPCA.mu for point in data_raw])                   
# rotate the centered data                                                      
centered_and_rotated = np.array([np.dot(myPCA.Wt, point) for point in centered])
# the variance along axis 0 should now be greater than along 1, so on           
variances = np.array([np.var(centered_and_rotated[:,i]) for i in range(3)])     
# they are not:                                                                 
print(variances[1]>variances[2]) #False; I want this to be True                 

# Now look at the PCA output, Y. This is centered, scaled, and rotated.          
# The variances decrease in magnitude, as I want them to:                       
variances2 = np.array([np.var(myPCA.Y[:,i]) for i in range(3)])                 
# This looks good, but the coordinates have been scaled.      
# Let's try to get from the raw coordinates to the PCA output Y
# mu is the vector of means of the raw data, and sigma is the vector of 
# standard deviations of the raw data along each coordinate direction             
guess = np.array([np.dot(myPCA.Wt, (xxx-myPCA.mu)/myPCA.sigma) for xxx in data_raw])
print(guess==myPCA.Y) # all true  

The last two lines in the above show that we may take a point $\mathbf{x}$ from its representation in the raw data input into its representation $\mathbf{x}'$ in terms of the PCA axes via $$ \mathbf{x}' = \mathrm{R}\cdot\left((\mathbf{x}-\boldsymbol{\mu}) / \boldsymbol{\sigma} \right) $$

where $\mathrm{R}$ is myPCA.Wt, the weight matrix, $\boldsymbol{\mu}$ is the vector of means of the original data along each coordinate axis, $\boldsymbol{\sigma}$ is the vector of standard deviations of the original data along each coordinate axis, and the division is element-wise. In order to write this in standard mathematical notation, let's replace this division by multiplication: $$ \mathbf{x}' = \mathrm{R}\cdot\left(\mathrm{D}\cdot(\mathbf{x}-\boldsymbol{\mu}) \right) $$ where $\mathrm{D}$ is a diagonal matrix whose diagonal entries are $1/\sigma_i$.

This notation makes clear the problem: to undo the scaling, I need to act on the RHS above with $\mathrm{R}\mathrm{D}^{-1}\mathrm{R}^{-1}$. This will return me to the problem situation, in which the variance is greater along the $z$ axis than the $y$.

Is there a way to use PCA to get what I want, or do I need to use another method?

  • $\begingroup$ I don't know about how Python does things, but doesn't the axis ordering depend on how NumPy's diagonalization procedure orders the singular values? Typically programs will either have them in ascending or in descending order of modulus. To that end, is there some reason you can't just reorder the axes to make $x$ and $y$ the largest variance? Or am I misunderstanding what you're asking? $\endgroup$ – DumpsterDoofus Aug 12 '14 at 20:32
  • $\begingroup$ @DumpsterDoofus The ordering is largest to smallest variance. This is the case for myPCA.Y. My question is why this is no longer the the case when I follow my procedure. I am reticent to reorder them because to do so requires reordering the columns of the rotation Wt, which could cause it to become a rotation and reflection rather than just a rotation (i.e., det(Wt)=1 and reordering can cause det(Wt)=-1). I could of course fix this by changing the signs on the columns, but this will result in reflections of the data, which I want to avoid $\endgroup$ – Eric Kightley Aug 12 '14 at 20:42
  • $\begingroup$ Ok, I think I sort of see what you're saying now. I think the problem is due to the fact that PCA is highly dependent on how your experimental variables are scaled. Quoting Wikipedia, "But if we multiply all values of the first variable by 100, then the first principal component will be almost the same as that variable, with a small contribution from the other variable, whereas the second component will be almost aligned with the second original variable." $\endgroup$ – DumpsterDoofus Aug 12 '14 at 22:27
  • $\begingroup$ More to the point, if you take the principle components computed by one application of PCA and then apply some sort of scaling to the original variables that make up those principal components, then (I think?) the new principal components are no longer uncorrelated, so you're sort of destroying the decorrelation that PCA is meant to give in the first place. Again quoting the Wiki, "This means that whenever the different variables have different units (like temperature and mass), PCA is a somewhat arbitrary method of analysis." So if scaling is arbitrary, PCA may not be the best tool. $\endgroup$ – DumpsterDoofus Aug 12 '14 at 22:36
  • $\begingroup$ Since $RD^{-1}R^{-1}$ is essentially just applying a scaling to the original variables that make up the principal components like I mentioned above, that's probably why you're running into the problem you're having. That said, I'm not experienced in statistics, so I'd ask others about what might be a useful workaround. $\endgroup$ – DumpsterDoofus Aug 12 '14 at 22:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.