Does Principal Component Regression still work in high-dimensional ($N<p$) situation?

I understand that, many classical methods for multiple regression won't work when $N<p$, where $p$ is the dimension of the input space and $N$ is the sample size.

For example, LSE for multiple regression, if $\mathbf{X}^T\mathbf{X}$ is nonsingular, then the unique solution is given by $\hat\beta=(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y}$. But if $p>N$, $\mathbf{X}_{N\times p}$ cannot be of full column rank, thus $\mathbf{X}^T\mathbf{X}$ is singular and the parameters in the regression model cannot be uniquely estimated.

Now, my confusion is, Principal Component Regression is also a popular classical regression method for high-dimensional data, but will there be any problem, like above, when we apply PCR to a high-dimensional ($n\ll p$) data?

More specifically, if $\mathbf{X}^T\mathbf{X}$ is singular, how to calculate the principal components? (How to make the eigen decomposition for a singular matrix?) Are the principal components given by softwares still correct? Can we still regress on these principal components for dimension reduction purpose?

Could anybody give me a hint? Thanks a lot!

-

With a rank-3 random correlationmatrix r $$\text{ R =} \small \left[ \begin{array} {rrrrrr} 1.0000& -0.8790& -0.8389& -0.8184& -0.5239& -0.7192\\ -0.8790& 1.0000& 0.6751& 0.5508& 0.0846& 0.3115\\ -0.8389& 0.6751& 1.0000& 0.9721& 0.7127& 0.7833\\ -0.8184& 0.5508& 0.9721& 1.0000& 0.8534& 0.9053\\ -0.5239& 0.0846& 0.7127& 0.8534& 1.0000& 0.9633\\ -0.7192& 0.3115& 0.7833& 0.9053& 0.9633& 1.0000 \end{array} \right]$$ I get by the "robust" cholesky-decomposition the loadings-matrix $L$ with only 3 significant columns $L = cholesky(R)$ $$\text{ L =} \small \left[ \begin{array} {rrrrrr} 1.0000& 0.0000& 0.0000& 0.0000& 0.0000& 0.0000\\ -0.8790& 0.4767& 0.0000& 0.0000& 0.0000& 0.0000\\ -0.8389& -0.1308& 0.5283& 0.0000& 0.0000& 0.0000\\ -0.8184& -0.3536& 0.4530& -0.0000& 0.0000& 0.0000\\ -0.5239& -0.7886& 0.3218& 0.0000& -0.0000& 0.0000\\ -0.7192& -0.6727& 0.1742& -0.0000& -0.0000& 0.0000 \end{array} \right]$$ and by iterative columnwise jacobi-rotation on the first 3 nonzero axes $P=\operatorname{rot}(L,\text{"PCA"})$ the principal components-solution $$\text{ P =} \small \left[ \begin{array} {rrr} 0.9037& -0.3761& 0.2045\\ -0.6478& 0.7610& -0.0363\\ -0.9557& 0.0870& 0.2812\\ -0.9832& -0.1061& 0.1482\\ -0.8119& -0.5821& -0.0446\\ -0.9087& -0.3731& -0.1871 \end{array} \right]$$ At the moment I do not exactly how to proceed to the regression-step - perhaps please give a hint how your model of dependent/independent variables is thought to be organized.
What I mean by Principal Component Regression is: suppose we have $\mathbf{X}^T\mathbf{X}=\mathbf{V} \mathbf{D}^2\mathbf{V}$, principal component regression forms the derived input columns $\mathbf{z}_m=\mathbf{X}v_m$ and then regresses $\mathbf{y}$ on $\mathbf{z}_1, \cdots, \mathbf{z}_M$ for some $M\leq p$. Because $\mathbf{z}_m$ are orthogonal, This regression is just a sum of univariate regressions. – Tansy Feb 2 '13 at 18:59