PCA when # of features larger than # of samples? [duplicate]

I'm using this code that adopts PCA for domain adaptation. I have some questions regarding how it handles the situation when the # of features is larger than the # of samples. Here is the snippet of the relevant Matlab code.

The input to the function is a design matrix $X$, where the rows correspond to observed samples and columns to features.

[nSmp,nFt] = size(X);

% centering X (zero-meaning)
mu = mean(X,1);
X = bsxfun(@minus,X,mu);

% select the faster way to do PCA
if nFt > nSmp
[V1, D1] = eig(X*X'/(nFt-1)); %%%%%%% Why?? %%%%%%
else
[V1, D1] = eig(X'*X/(nFt-1)); %% standard doing %%
end

% sort the eigenvalues and re-order the eigenvectors
% in descending eigenvalues order
D1 = diag(D1);
[evs,I] = sort(D1,'descend');
V1 = V1(:,I);

% cumulative percentage of variance accounted
D1 = cumsum(evs);
D1 = D1/D1(end);

% W_prj is the pca coefficient
if nFt > nSmp
W_prj = X'*V1; %%%%%% Why?? %%%%%%%
else
W_prj = V1;    %% standard doing %%
end

% normalize PCA coefficient so that for each column, the L2 norm is 1
W_prj = bsxfun(@rdivide,W_prj,sqrt(sum(W_prj.^2,1)));

% project X to PCA subspace
Xnew = X*W_prj;


I do not understand the parts where I put a "why??" in comment.

I know that eigendecomposing $X^TX$ is the standard way of computing PCA coefficients, as is done in the else clause of the first if statement. The PCA coeffcients are eigenvectors that are reordered according to the eigenvalues.

But when the # of feature is larger than # of samples, why is it OK to first eigendecompose $XX^T$, reorder the eigenvectors according to eigenvalues, and then return $X^TV_1$ as the PCA coefficients? Are there some kernel methods involved in this doing? What's the rationale behind this doing?

I have asked the same question here. When either post gets an answer, the other post will be deleted.

It's because: For an eigenvalue $\lambda$ and an eigenvector $\textbf{v}$, if $XX^T\textbf{v}=\lambda\textbf{v}$, then $X^TXX^T\textbf{v}=\lambda X^T\textbf{v}$, i.e. $X^T\textbf{v}$ is the eigenvector of $X^TX$.
As we all know that eigendecomposition is a time-consuming process, by this way we can always decompose a small matrix whose size is $\min\{n\_feature, n\_sample\}$.