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?
