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?
 A: Here is a relevant reference:
How to generate random matrices from the classical compact groups
http://arxiv.org/abs/arXiv:math-ph/0609050
A: 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.
A: The best algorithm is the one due to G.W. Stewart
Stewart, G. W., The efficient generation of random orthogonal matrices with an application to condition estimators. (With mircofiche section), SIAM J. Numer. Anal. 17, 403-409 (1980). ZBL0443.65027.
Algorithm: Let $M$ be $n\times n$ with i.i.d $N(0, 1)$ entries. Let $M=QR$ be the QR factorization ($Q$ orthogonal, $R$ upper triangular). Then $Q$ is uniformly distributed on $O(n)$ (project to $SO(n)$ if you like).
