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

## 1 Answer

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