# 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.

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$)

• Wow this was beautiful. Thank you for this excellent explanation of precisely what was going on. – The_Anomaly Aug 16 '17 at 19:48