Sample a random rotation in $n$ dimensions 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.
 A: 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

A: 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)$.
A: 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.
