Due to the idempotence property of projections that you quote, a projection is its own unique pseudoinverse, as may be shown by substituting in conditions (I) to (IV) of the definition of a Moore-Penrose pseudoinverse. Other generalised inverses satisfying properties (I) and (II) will in general exist, but the MP pseudoinverse is optimal in a minimum norm sense, that is to say there is nothing better. You already have the least squares solution you are looking for I'm afraid.
Multiplying a matrix by a generalised inverse will not in general give the identity matrix on either side, viz. it may be lacking in rank on both sides, and this is in general the case for projections. The identity matrix is the exception as you observe.
You are correct in suggesting that the lost dimensions cannot be recovered.
I wouldn't say we "prefer" the eigenvalues to the singular values in the symmetric/Hermitian case. They are effectively the same thing! Or, more accurately, the singular values are nothing more than the absolute values of the eigenvalues.
To see why, let $X$ be any Hermitian matrix. Its Schur (eigenvalue) decomposition is $X=Q\Lambda Q^H$ where $\Lambda$ is diagonal, and $Q^HQ=QQ^H=I$. Its singular value decomposition, on the other hand, is $X=U\Sigma V^H$, where $\Sigma$ is diagonal and nonnegative, and $U^HU=UU^H=I$ and $V^HV=VV^H=I$.
Alternatively, we can write these decompositions in dyadic form:
$$X=\sum_{i=1}^n \lambda_i q_i q_i^H = \sum_{i=1}^n \sigma_i u_i v_i^H$$
where $q_i$, $u_i$, and $v_i$ are the $i$th columns of $Q$, $U$, and $V$, respectively, and $\lambda_i$ and $\sigma_i$ are the $i$th diagonal elements of $\Lambda$ and $\Sigma$, respectively.
In this latter form, it is simple to see that, given the Schur decomposition,
$$\sigma_i = |\lambda_i|, \quad u_i = q_i, \quad v_i = \begin{cases} v_i & \lambda_i \geq 0 \\ -v_i & \lambda_i < 0 \end{cases}
\qquad i=1,2,\dots, n$$
gives us a valid singular value decomposition. Our choice to apply the sign flipping to $v_i$ is arbitrary.
Be assured that in practice, this is precisely what you should be doing to compute the singular value decomposition of a Hermitian matrix. The symmetric eigenvalue problem is well-conditioned, and there is no numerical advantage to applying the a nonsymmetric SVD to a Hermitian matrix. In fact, the computation of the singular values of a non-Hermitian matrix $A$ are sometimes computed by applying a symmetric eigenvalue algorithm to the Hermitian matrix
$$\begin{bmatrix} 0 & A \\ A^T & 0 \end{bmatrix}$$
I'll leave it as an exercise for you to prove that this works :-)
EDIT: DGrady's comment alerted me to the fact that I didn't read the second half of the question correctly. In it, he was comparing the merits of computing the SVD of $X$ versus the eigenvalue decomposition of $C=XX^H$. In this case, it is indeed going to be a bit better to compute the SVD of $X$ directly.
Best Answer
$M^+=V\Sigma^+U^\ast$ (i.e. $M^+=V\Sigma^+U^\top$ in the real case), where $\Sigma^+$ is a rectangular diagonal matrix whose size is identical to the size of $\Sigma^\top$. The $i$-th main diagonal entry of $\Sigma^+$ is $\sigma_i^{-1}$ if the $i$-th singular value $\sigma_i$ of $M$ is positive, otherwise the diagonal entry is zero.
You may simply prove that $M^+$ is indeed the Moore-Penrose pseudoinverse of $M$ by showing that it satisfies the four defining properties of Moore-Penrose pseudoinverse, namely, both $MM^+$ and $M^+M$ are Hermitian (or real symmetric in your case), $MM^+M=M$ and $M^+MM^+=M^+$.