# How to correctly generate uniformly distibuted random elements from SO(n)?

I already found some way to produce such matrices from SO(n) with a method called subgroup algorithm but I would like some advice on the method I used. Nowhere I could really find any paper relating to SO(n) directly but rather O(n).

To make it simple, I generate a $(k-1) \times (k-1)$ rotation matrix $\Gamma_{k-1}$, create the matrix

$\left( \begin{array}{ccc} 1 & \cdots & 0 \\\ \vdots & \Gamma_{n-1} & \\\ 0 & & \end{array} \right)$

then I find a rotation $R$ that transforms the first column into a vector $v$ generated randomly as a point on the unit $k$-sphere. Then generate a new $k \times k$ rotation matrix $\Gamma_k$ :

$\Gamma_{k} = R\left( \begin{array}{ccc} 1 & \cdots & 0 \\\ \vdots & \Gamma_{k-1} & \\\ 0 & & \end{array} \right)$

Induction is doing all the work. In practice, either you start by a $1 \times 1$ matrix with single element equal to 1, or you create a first trivial $2 \times 2$ rotation matrix with a angle taken uniformly from $[0,2\pi)$ which seems more direct. I skipped the details on how to generate the matrix $R$ and the vector $v$ as it seemed useless at first. To sum up $v$ is the normalization of a vector $v'$ which elements $v'_i \sim N(0,1)$ and $R$ is obtained by a double Householder reflection.

In practice this seems like a correct incarnation of the subgroup algorithm where the group is SO(n) but I wanted to check with more mathematically inclined than me if it didn't exist a more efficient way! As I said, I found some papers explaining variations or improvements on the original method from Stewart$-$which the subgroup algorithm is a generalization$-$in papers from Anderson or Genz. All of them focused on orthogonal matrices and I would have liked to transpose these improvements to special orthogonal matrices.

See:

• "The Efficient generation of Random Orthogonal Matrices with an Application to Condition Estimators" Stewart
• "Generation of Random Orthogonal Matrices" Anderson
• "Methods for Generating Random Orthogonal Matrices" Genz
• "The Subgroup Algorithm for Generating Uniform Random Variables" Diaconis

EDIT: by the way, I might use this algorithm with n as large as 1000 and beyond

So the real question is (and I'm not really sure how to formulate that correctly as I have a serious lack of knowledge in that area), is there a mapping from O(n) to SO(n) that will preserve the uniformity?

-
If the resulting matrix is in $O(n)$ but not in $SO(n)$ simply exchange two columns. An orthogonal matrix is one whose columns define an orthonormal basis. The difference between $O(n)$ and $SO(n)$ is just orientation, so reorder the basis to get the right orientation. This of course is independent of the algorithm you choose to generate the element. – José Figueroa-O'Farrill Oct 13 '11 at 15:22
Is there something about this algorithm you don't like? Choose all entries of a $n\times n$ matrix with a pseudorandom number generator. Then apply gram-schmidt. – Ryan Budney Oct 13 '11 at 19:24
@Ryan: Gram-Schmidt does not yield a uniform distribution. – S. Carnahan Oct 13 '11 at 20:14
@Scott: which version are you using? There are symmetric versions of Gram-Schmidt. – Ryan Budney Oct 13 '11 at 20:23
If there is a version of Gram-Schmidt that yields a uniform distribution, then I think the choice of version is rather important, and should be precisely specified. Even so, the runtime may have unfavorable asymptotics (but I'm on shaky ground here). – S. Carnahan Oct 13 '11 at 22:00

Here is a relevant reference: How to generate random matrices from the classical compact groups http://arxiv.org/abs/arXiv:math-ph/0609050

-
Nice reference... – Igor Rivin Oct 13 '11 at 19:21

You can extract an answer from the Wikipedia article (and the Diaconis-Shahshahani paper that is referenced). First, find/write a function that yields a uniform distribution of points on a unit sphere (e.g., divide a normal distribution in Euclidean space by norm). Then, if you have a uniform distribution on $O(n-1)$, embedded in $O(n)$ in the way you described, you can just reflect in (the hyperplane orthogonal to) your random unit vector to get a uniform distribution on $O(n)$.

If you want a uniformly distributed element of $SO(n)$, you should start your process using the 1-by-1 matrix $(-1)^{n-1} \in O(1)$ instead of a random element of $\{ -1, 1\}$. That will fix the determinant.

Edit: Diaconis-Shahshahani prove that the subgroup method that I described above yields a distribution that is uniform with respect to the Haar measure on $O(n)$. If you are getting strange things from a visual check, then the problem is almost certainly with your implementation.

The method I give for obtaining a uniform distribution on $SO(n)$ is a straightforward alteration of the Diaconis-Shahshahani method. At each step of their recursion, the determinant of your matrix is multiplied by $-1$. Therefore, the determinant you get at the end depends only on the dimension of the final matrix, and the determinant of the initial $1 \times 1$ matrix. If you choose this initial matrix to be $1$ or $-1$ and iterate the process as before, this will concentrate the distribution on to one of the cosets of $SO(n)$ in $O(n)$, but it will be uniform with respect to the translation action of $SO(n)$. Therefore, you only need to choose the correct initial matrix (using the formula I provided) to get a uniform distribution on $SO(n)$ by this method.

An alternative method was given by José in the comments. Composition with any determinant $-1$ element takes the uniform distribution on $SO(n)$ to the uniform distribution on the nontrivial coset in $O(n)$, and back. Do the subgroup method recursion with random initial step, and if at the end the determinant is $-1$ (which will happen precisely when your dimension is even), switch a pair of rows. Any pair of rows will do.

Here is the answer to your latest question: the map that switches the first two rows of a matrix if and only if the determinant of that matrix is $-1$ is a mapping from $O(n)$ to $SO(n)$ that preserves uniformity.

-
I used that paper but again it talk about O(n) but not SO(n), I didn't find the method you give in this form. I doubt the resulting matrices are distributed with Haar measure. However, I barely scratched the surface of that recently and delved into literature... my curriculum was lacking some advanced math courses. In 3d this seems to work well and might even be correctly distributed, but as 3d is a particular case... I checked visually and I saw strange things going on! – LeftBrainDamaged Oct 14 '11 at 14:21
Many thanks that will certainly help a lot! Would it be possible that the problem doesn't lie in my implementation but rather in the choice of hyperplane for reflections? You use directly the random unit vector as the normal for this hyperplane, which seems causing the strange visual behavior in my tests... – LeftBrainDamaged Oct 25 '11 at 10:00
A uniform distribution of unit vectors gives you a uniform distribution of Householder reflections. – S. Carnahan Oct 26 '11 at 2:37
S. Carnahan answer is actually wrong because of one thing: the Diaconis-Shahshahani paper redefine the Householder reflection to make it unique so as to get an element of O(n) distributed with Haar measure (to get uniformity). Classical Householder transformation is a reflection and thus determinant is toggled between -1 and 1 at each iteration of Diaconis' algorithm, however this is not always the case with the alteration. As a result you cannot rely anymore on this toggling and use the first algorithm S. Carnahan described. The alternative method -- switching a pair of rows (or columns) of t – LeftBrainDamaged Apr 9 '12 at 16:32