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?**