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.