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?


$P$ won't necessarily be a rotation matrix.

To see why, let $p = 3$, suppose we generate the lower-triangular part of the matrix $P$ according to this scheme as $$ P = \left[ \begin{matrix} * & * & * \\ p_{21} & * & * \\ p_{31} & p_{32} & * \end{matrix} \right] $$ There is no guarantee that $p_{21}^2 + p_{31}^2 < 1$, and hence no guarantee that there will exist a value of $p_{11}$ that will make the first column a unit vector, because $p_{11}^2 + p_{21}^2 + p_{31}^2 \ge p_{21}^2 + p_{31}^2 $. In fact, simple geometry tells us that $p_{21}^2 + p_{31}^2 < 1$ will happen only with probability $\frac{\pi}{4}$ (think of the unit square and unit circle overlayed), and thus there is a $1 - \frac{\pi}{4} \approx 21\%$ chance that this algorithm will fail in this case.

What the authors might have meant was to generate $P$ by drawing from $U(0, 1)$ and then by applying the Gram-Schmidt process, as described in this math.stackexchange question, but this is still the "wrong" way to generate random orthogonal matrices if what you really want are uniformly random orthogonal matrices, as is described in the accepted answer to that question.

| cite | improve this answer | |
  • $\begingroup$ Thanks, I was hoping for an easier way, but that's life :) $\endgroup$ – mariob6 Sep 8 at 6:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.