Let us consider the Riemannian manifold $$M = \{ O \in \mathbb R^{n \times n}\mid O^T O = I \}$$ with inner product $$\langle P, Q \rangle_M \equiv \mbox{Tr} \left( P^T Q \right)$$ which naturally embeds into the ambient space $\mathbb R^{n \times n}$ via the identity embedding. It can be shown that the tangent space at a point $P \in M$ is $$T_P M \equiv \{ PZ \mid Z \in \mathbb R^{n \times n}, ~~ Z^T + Z = 0 \}$$
Given a general matrix $N \in \mathbb R^{n \times n}$, how does one compute the projection of $N$ onto $T_p M$, with $T_p M$ interpreted as a subspace embedded in $\mathbb R^{n \times n}$?
In the paper The geometry of algorithms with orthogonality constraints, this is claimed to be:
$$\Pi_{T_PM}(X) = P (P^T X - X^T P)/2 + (I - PP^T)X \tag{eqn 2.4}$$
I can verify that this is correct (as is done in this math.SE answer). The question is how to derive this for the $O(n)$ case, in a way that generalizes to other matrix manifolds?
My attempt at an answer
We can try to take an orthornormal basis of skew-symmetric matrices $\{ K_i \}$, and then use these to create a basis for $T_p M$ as $\{ P K_i \}$.
So this leads to the ugly formula:
\begin{align*} \Pi_{T_PM}(X) \equiv \sum_i \mbox{Tr}(X^T P K_i) P K_i = P \sum_i \mbox{Tr}(X^T P K_i) K_i \end{align*}
I'm unable to simplify the expression any further, unfortunately. I presume I am missing some piece of structure of the skew-symmetric matrices, or some tool from manifold/lie theory to close the proof.
I will assume you are interested in the projection onto matrix subgroups. Its really only necessary to find this projection at the identity. If $\pi: T_1\Bbb R^{n\times n}\to T_1 M$ is the orthogonal projection at $1$, then $P\cdot \pi\cdot P^{-1} : T_P\Bbb R^{n\times n}\to T_PM$ is the projection at $P$.
In your case it is very natural that the map $X\mapsto \frac{X-X^t}2$ is the projection onto the anti-symmetric matrices (ie to $T_1O(n)$), hence at a general $P$ you recover $X\mapsto \frac{X-PX^tP}{2}$, which agrees with your formula.
More generally suppose that your matrix group is of the form $G=\{ A\mid A^tJA=J\}$ for some invertible matrix $J$ (this covers a very large amount of matrix groups). Then you may calculate that the tanget space at $1$ is $T_1G= \{ X\mid X^tJ+JX=0\}$. The natural formula for a projection to this space is then $X\mapsto \frac{X-J^{-1}X^t J}2$. I'll leave the verification that this is an orthogonal projection to you (I haven't checked).
Additionally lets look at the case $G=SL_n$. Here the group condition is $\det(A)=1$, which translates to $\mathrm{Tr}(X)=0$ for the Lie algebra. The map $X\mapsto X-\frac1n\mathrm{Tr}(X)$ is the natural candidate for the projection here.