5

How do I sample a random rotation matrix in $n$ dimensions?

Let $M$ be a $n\times n$-dimensional matrix. $M$ is a rotation matrix if $\|Mx\|_2 = \|x\|_2$ for all $x \in \mathbb{R}^n$. I want to sample uniformly at random from the rotations, in the following sense. Let $\mathcal{S}=\{x: \|x\|_2=1\}$ denote the perimeter of a sphere of radius $1.$ I want a distribution on $M$ so that, if you fix any point $x \in \mathcal{S}$, then $Mx$ will be uniformly distributed on $\mathcal{S}$.

How can I sample a matrix $M$ from this distribution? I am looking for an algorithmic procedure I can use to sample such a matrix randomly.

I believe this corresponds to sampling uniformly at random from $O(n)$ (the orthogonal group). I am fine with either a procedure for sampling uniformly at random from $O(n)$ or $SO(n)$; it should be easy to convert between the two of them.

Michael Hardy
  • 1
D.W.
  • 5,958
  • Your definition also allows reflections (elements of the orthogonal group $O(n)$ of determinant $-1$) but I suppose the difference isn't too bad either way. – Qiaochu Yuan Sep 24 '20 at 20:02
  • @QiaochuYuan, oh good point. I'm happy with either a procedure for sampling from $O(n)$ or $SO(n)$. Thank you for your answer! – D.W. Sep 24 '20 at 20:05
  • 1
    In the case of $n=3$, a uniformly distributed rotation can be derived efficiently from a uniformly distributed point on the 3-sphere $S^3\subset\mathbb{R}^4$ via a simple formula. – Kajelad Sep 24 '20 at 23:59
  • 1
    See here for the case $n=3$, and a handful of ideas... – Jean Marie Jul 23 '24 at 22:03

3 Answers3

4

Here's a fun indirect way to do it: sample a random $n \times n$ real matrix $M$ from the Gaussian orthogonal ensemble. Explicitly this means that the diagonals $M_{ii}$ are independent Gaussians $N(0, \sqrt{2})$ and the off-diagonals $M_{ij} = M_{ji}$ are independent Gaussians $N(0, 1)$ subject to the constraint that $M$ is symmetric. Now compute the eigendecomposition $M = UDU^{\dagger}$. Then $U \in O(n)$ will be distributed according to Haar measure. (Maybe it would be better to use the QR decomposition, I don't know.)

There are efficient methods for both generating a bunch of independent Gaussians and for computing eigendecompositions in, say, Octave or MATLAB so this can be done with very little code, although I don't know what the asymptotics are or how much better it's possible to do. In Octave it is something like this (untested):

A = randn(n); # generates an nxn matrix all of whose entries are N(0, 1)
M = (A + A.') / sqrt(2); # symmetrizes A, now a sample from GOE
[V, E] = eig(M); # computes eigenvectors V and eigenvalues E of M
Qiaochu Yuan
  • 468,795
  • Hm, combining this recipe with "swap first two rows to make determinant positive", distribution of angles appears unsatisfyingly non-uniform ... https://math.stackexchange.com/questions/4962936/sampling-from-so4-more-uniformly-than-haar – Yaroslav Bulatov Aug 25 '24 at 16:58
  • @Yaroslav: then don't swap the first two rows to make the determinant positive! – Qiaochu Yuan Aug 25 '24 at 19:03
4

As a more direct (and, I suspect, much more computationally efficient) version of Qiaochu Yuan's approach, we can bypass the need for a diagonalization with a slightly simpler ensemble, and use the fact that the Graham-Schmidt process is $O(n)$-invariant.

Define a random $n\times n$ matrix $M$ where each entry is an iid. Gaussian. This ensemble is $O(n)$ invariant, in the sense that $p_M(A)=p_M(OA)$, where $p_M$ is the probability distribution of $M$ and $O\in O(n)$. To avoid numerical issues, you can freely reject matrices with $|\det(M)|<\epsilon$ for some fixed $\epsilon$, and still obtain an $O(n)$-invariant ensemble since the absolute value of the determinant is $O(n)$-invariant.

Applying the Graham-Shmidt process $\Gamma$ to $M$ yields a new random matrix $\Gamma(M)\in O(n)$. Since $\Gamma(OA)=O\Gamma(A)$ for $O\in O(n)$, $\Gamma(M)$ is also $O(n)$ invariant, which is to say it is uniformly distributed w.r.t. the Haar measure on $O(n)$.

Kajelad
  • 15,356
  • 1
    $SO(n)$ was a typo. This process doesn't control the sign of the determinant, so it would be uniform on $O(n)$. It's not difficult to restrict to a uniform sampling of $SO(n)$, though, e.g. by flipping the sign of the first column if needed. – Kajelad Sep 26 '20 at 21:27
3

According to a similar question in stack overflow, in Python there is special_ortho_group in scipy.stats (documentation), which you can use as follows (copied from stack overflow)

from scipy.stats import special_ortho_group
num_dim = 3
x = special_ortho_group.rvs(num_dim)

I'm adding this answer in case someone is searching for a simple computer implementation rather than a mathematical understanding of it. For the mathematics, either check the other answers or the source code.

Carlos Pinzón
  • 1,333