0
$\begingroup$

I am learning PCA, and while this seems like an obvious question, I can't seem to understand some of the ideas behind which eigenvectors to use when performing orthogonal regression with PCA. My code is below, and a diagram.

I am using a parametric model of M=[X.^2;X;Y] to estimate a quadratic with parameters [2;0.7;2]. I am performing eigenvalue decomposition on that. In doing so, I am finding out the directions of variance in the data, and the scale of such variance. PCA aims to find which features best describes the data?

So when using the below code, I find that estimating with the first two eigenvectors (directions) (V(:,1:2)') provides a good fit of the function (shown by the green circles in the diagram). However, I cannot get a parameter estimate ([2;0.7;2]) from this data (maybe you can, but I couldn't figure out how). So instead I looked at the final eigenvector (V(:,3)), and this gives me the parameters of my function (by taking V(1:2,3)/V(3,3)). Applying this results in an estimation as shown by the red crosses.

As you can see from the diagram, both visual fits are exact. My question is, how can I have estimated this function in two separate ways? I thought that if I used an eigenvector that did not have the smallest eigenvalue, then I would be finding a line in some other direction, not the main direction?

Any help anyone can give me would be appreciated, I think I understand how to apply PCA, just unsure as to this directionality thing.

% PCA on Quadratic using parametric model

% Setup data
X=linspace(-1,1,100);
Y=2*X.^2+0.7*X+2;
M=[X.^2;X;Y]';

% Normalize around centroid
meanM=mean(M);
MAdjust = M-repmat(meanM,100,1);

% Perform eigenvalue decomposition
[V D]= eigs(cov(MAdjust));

% Estimate function using V(:,1:2)
FinalData1 = V(:,1:2)'*MAdjust';
MHat1=repmat(meanM,100,1)+(pinv( V(:,1:2))'*FinalData1)';

% Plot starting data, and prediction using V(:,1:2)
figure;hold on;
scatter3(M(:,1),M(:,2),M(:,3),15,'k');
scatter3(MHat1(:,1),MHat1(:,2),MHat1(:,3),50, 'g','d');

% Estimation function using V(:,3)
c = meanM*(V(:,3)/V(3,3));
param = -V(1:2,3)/V(3,3)
yhat =  [X.^2;X;X.^0]'*[param ; c];

% Plot the estimated function
scatter3(X.^2,X,yhat,50, 'r','+');

xlabel('X^2');ylabel('X');zlabel('Y');

3D plot of 2*X.^2+0.7*X+2 with predictions from eigenvalue decomposition parameter estimates

$\endgroup$
1
$\begingroup$

I am coming late at the party but do you have a particular reference you used for your 3rd step (Estimation function using V(:,3))? To refer to your question exactly, you did not "estimate this function in two separate ways". In the first case you took the projection of your sample using the two PCs that fully encapsulate all the variance in your sample. No estimation took place. The omitted eigenvector had an associated eigenvalue of $0$. In the second step (the one I am not actually clear where you are based) you estimated the last column of your data as a linear combination of your previous columns.

Some general comments: 1. Based on the eigenvalues you receive, this is a 2-D dataset so no dimension reduction takes place. 2. The transpose of the pseudo-inverse of orthogonal vectors (pinv( V(:,1:2))') are just the vectors you started with (V(:,1:2)); you simply lose precision with this calculation. 3. Use princomp to simply your life a bit. 4. You are estimating the covariance matrix but also include $Y$. That seems a bit odd.

I believe what you want, given $X^T X = V D V^{T}$, is : $W = XV$ and then $\beta^* = (W^T W)^{-1} W^T Y$ and $\beta = V \beta^*$. This will get you the parameter estimates you inquire for. So for example:

x=linspace(-1,1,100);
Y=2*x.^2+0.7*x+2;
X = [ x'.^2  x' ones(100,1)]; 
[V,D] = eigs(cov(X)); % Step 1
W = X * V'; % Step 2
beta_PCR = W \ Y' % Step 2 cont. 
% beta_PCR =
%   0.700000000000000
%   2.000000000000000
%   2.000000000000000
beta_original = V * beta_PCR % Step 3 
% beta_original =
%   2.000000000000000
%   0.700000000000000
%   2.000000000000000

yhat_1 = beta_original' * [ x'.^2   x' ones(100,1)  ]';

In broad strokes the principal component regression procedure described above does the following steps:

  1. PCA on the matrix $X$ where the columns of $X$ are the independent variables (getting relevant PCs).
  2. Regression between $Y$ and independent variables projected in the subspace defined by the PCs selected (getting relevant $\beta$'s on that domain).
  3. Projection of the $\beta$'s calculated in the previous step back to their original scale using the PCs selected in at the beginning of this procedure.

Note: You might be tempted to use two eigenvectors in the example above based on the fact that one eigenvalue is $0$. This is not the same as with the case where $Y$ was included. The $0$-th eigenvalue here is because we include a vector of $1$'s for the intercept; clearly a stable vector is not covarying with anything. Additionally, the data "are centred" when we do the PCA as we do use the covariance matrix. We do not centre $X$ when projecting them as we want to keep $\beta$ interpretable. If we centred $X$ during the project we would need to have something like:

newX = X - repmat(mean(X),100,1);
[V,D] = eigs(cov(newX));
W = newX * V';
beta_PCR = W \ Y'
%beta_PCR =
%  -0.699999999999999
%   1.999999999999997
%                   0
beta_original = V * beta_PCR
%beta_original =
%  -1.999999999999997
%  -0.699999999999999
%                   0

yhat_2 = beta_original' * newX' - mean(Y);

which while equally correct with the previous expression does not allow the direct interpretation of the coefficient you get. Plotting the estimated values confirms the equivalence of the two approach easily:

plot(Y,'b-','LineWidth',5); hold on; grid on;
plot(-yhat_2,'r--', 'LineWidth',4); plot(yhat_1,'ks')
legend('Original data', 'Est. Values #1', 'Est. Values #2')

enter image description here

$\endgroup$

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.