I found the recipe below to work for generating random rotation matrices in n dimensions. (straightforward modification from CUE sampling recipe from page 11 of https://arxiv.org/abs/math-ph/0609050). Is there a similarly efficient algorithm to generate small rotations? (ie, I want the vectors before and after transformation to be close in some sense)
randomRotation[n_]:=Module[{M,z,q,r,d,ph},
z=RandomVariate[NormalDistribution[0,1],{n,n}];
{q,r}=QRDecomposition[z];
d=Diagonal[r];
ph=d/Abs[d];
M=q*ph;
(* Note, determinant may be -1 giving "pseudo-rotation", switch 2 rows of the matrix to guarantee true rotation *)
indices=If[Det[M]>0,Range[n],{2,1}~Join~Range[3,n]];
M[[indices]]
];
You obtain Haar-uniform random rotations if you perform QR decomposition on an initial matrix which has independent and identically gaussian distributed random elements.
If you want random small rotations from QR, you could perform it on initial matrices which consist of random small perturbations of the identity.