I'm reading the article "Multivariate mixtures of normals with unknown number of components" (Dellaportas and Papageorgiou 2006).
In this article they describe in great details how to implement a Reversible Jump MCMC algorithm, similar to the one by Richardson and Green, when the data is living in $R^d$.
The key point is that one needs to randomly generate a $(d \times d)$ rotation matrix $P$ in order to propose new values for the split move. When $d=2$, this simply amounts to generate $\theta \in [0, \pi]$.
However when $d \geq 3$, the authors say:
Let also $P$ be a $(p × p)$ rotation matrix with columns orthonormal unit vectors which has $p(p − 1)/2$ free parameters. We generate $P$ by generating its lower triangular matrix under the diagonal independently from $p(p − 1)/2$ uniform $U(0, 1)$ densities
However, how can this $P$ be a rotation matrix?