Principal Component Regression: to Center or not to Center I am learning about Principal Component Regression, but I am confused about centering the matrix.  My question is described using the following example:
In MATLAB, I made some sample data, where $Y$ is a linear combination of the columns of $X$:
clc;clear all;
X = rand(1E3,20);
beta = rand(20,1);
Y = X*beta;

Of course, linear regression works perfectly:
%Compute the least squares using SVD
[U1, S1, V1 ] = svd(X,0)
beta0 = V1*inv(S1)*U1'*Y;
mean((X*beta0-Y).^2) %returns ~1E-30

Now, I will do PCR.  First, let's get the principal components of X, using SVD:
X_c = bsxfun(@minus, X, mean(X,1));
[U2,S2,V2] = svd(X_c);

%Principal components of X
W= X_c *V2;

And now, I get the regression coefficients of the principal components:
%Compute the regression coefficients of the PCs.  
[U3, S3,V3] = svd(W,0);
beta1 = V3*inv(S3)*U3'*Y;

But this linear combination of principal components does not recover Y.
Y2 = W*beta1;
mean((Y2-Y).^2) %returns ~20.1

However, if I use the uncentered matrix on the loadings of the principal components, it works nicely:
Y3 = X*V2*beta1; 
mean((Y3-Y).^2) %returns 1E-26

I am confused about how centering plays into this.  We compute the PCs (from the centered data), regress Y to the PCs, but then we need to project the original (uncentered data) using the loadings of the PCs, in order to use these regression coefficients.  Why does this work?  I would have thought that since the PCs were computed from the centered data, the regression coefficients would not be applicable to the uncentered data.
 A: The way you tried to perform principal component regression didn't work because you didn't include a constant term in the model. You defined $y = X \beta$, so regressing $y$ against $X$ without a constant term works in this case (but isn't a good idea in general; it won't work if $y = X \beta + c$ with nonzero $c$). You then computed the PCA projections $W$, which are a translation of $X$ (subtracting the mean), followed by a rotation. Because of the translation, $y$ can't be expressed as a linear combination of the columns of $W$. It will work if you include a constant term in the model: $y \approx W \hat{\beta} + c$
Why the last step worked
In the last step, you performed:
$$y_3 = X V_2 \beta_1$$
The way you computed $\beta_1$ using the SVD is equivalent to: $\beta_1 = W^+ y$, where $W^+$ is the Moore-Penrose pseudoinverse of $W$. Plug this in:
$$X V_2 W^+ y$$
You computed $W$ as: $W = X_c V_2$. Plug this in:
$$X V_2 (X_c V_2)^+ y$$
Rewrite this as:
$$X V_2 V_2^+ X_c^+ y$$
$V_2$ is orthonormal, so $V_2 V_2^+$ is the identity matrix:
$$X X_c^+ y$$
You defined $y$ as $y = X \beta$. Plug this in:
$$X X_c^+ X \beta$$
$X_c$ is a centered version of $X$, so we can write $X = X_c + \vec{1} \mu^T$, where $\vec{1}$ is a column vector containing $n$ ones (one for each data point), and $\mu^T$ is a row vector containing the mean over rows of $X$. Plug this in for the second instance of $X$:
$$X X_c^+ (X_c + \vec{1} \mu^T) \beta$$
Expand:
$$X X_c^+ X_c \beta + X X_c^+ \vec{1} \mu^T \beta$$
By the definition of the pseudoinverse, $X_c^+ X_c$ is the identity matrix:
$$X \beta + X X_c^+ \vec{1} \mu^T \beta$$
$X_c^+ \vec{1}$ can be re-written as $(\vec{1}^+ X_c)^+$. Note that $\vec{1}^+ X_c$ is the mean of $X_c$, which is zero. Therefore, the second term above reduces to zero and we're left with:
$$X \beta$$
which is equal to $y$ (in this case; as above, this won't work if $y = X \beta + c$)
