How to randomly sample small rotations?

Question

I need a random rotation matrix that rotates random vectors by $\approx \arccos 0.9$

• RandomVariate[CircularRealMatrixDistribution[n]] gives rotation matrices that are unrestricted
• Orthogonalizing IdentityMatrix[n] + small_noise gives a small rotation, but for larger $n$, this gives a large rotation even for small amount of noise.

Any tips?

Below is a heuristic attempt which adjusts amount of noise until arccos[0.9] target is hit, but it's too slow to run for $n=1000$

ClearAll["Global`*"];
randn[i_] := RandomVariate[NormalDistribution[], i];
randn[i_, j_] := RandomVariate[NormalDistribution[], {i, j}];

(*Generate rotation matrix by orthogonalizing Identity with eps \
perturbation*)
epsRotation[n_, eps_] := 
  Module[{M, z, q, r, d, ph, indices}, 
   z = IdentityMatrix[n] + eps   randn[n, n];
   {q, r} = QRDecomposition[z];
   d = Diagonal[r];
   ph = d/Abs[d];
   M = q*ph;
   (*determinant may be -1 corresponding to reflection,
   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]]
   ];

(* Finds eps such that average rotation is about arccos(0.9) *)
On[Assert];
getEps[d_] := Module[{},
   x0 = {1}~Join~ConstantArray[0, d - 1];
   target = 0.9;
   {epsHist, cosHist} = {{}, {}};
   eps = 0.001;
   For[i = 1, i < 1000, i += 1,
    cos = First[x0 . epsRotation[d, eps]];
    If[cos < target, eps *= 0.9, eps *= 1.1];
    AppendTo[epsHist, eps];
    AppendTo[cosHist, cos];];
   Assert[0.5 < Mean[cosHist[[-500 ;;]]] < 0.95];
   Mean[epsHist[[-500 ;;]]]
   ];

With[{d = 100},
 mat = epsRotation[d, getEps[d]];
 vec = randn[d];
 angle = vec . mat . vec/Norm[vec]^2;
 Print["rotation angle cos: ", angle]
 ]

Answer 1

The solution is sample a random element of the space of isoclinic rotations. Use the fact that real Schur form of orthogonal matrix is block diagonal with $2\times 2$ blocks corresponding to planar rotations. (8.2 of Handbook of Linear Algebra)

Sample eigenvectors from $O(n)$ (CircularRealMatrixDistribution), and use RotationMatrix[angle] for the $2\times 2$ blocks.

For $n\times n$ matrices, this space has $\frac{n(n-2)}{4}$ dimensions (proof)

randomIsoclinic[n_?EvenQ, theta_] := 
  Module[{blocks, realSchur, p, R},
   blocks = Table[R[i] -> RotationMatrix[theta], {i, n/2}]; 
   realSchur = ArrayFlatten[DiagonalMatrix[Array[R, n/2]] /. blocks]; 
   p = RandomVariate@CircularRealMatrixDistribution@n;
   p . realSchur . p\[Transpose]
   ];
angle = 1;
n = 4;
numSamples = 10;
mat = randomIsoclinic[n, angle];
data = Normalize /@ 
   RandomVariate[NormalDistribution[], {numSamples, n}];
AllTrue[data, ArcCos[# . mat . #] == angle &]

Answer 2

Could you not use spherical linear interpolation here in conjunction with RandomPoint[Sphere[n]]?

slerpAngle[v1_, v2_, a_] := 
 With[{Ω = ArcCos[v1 . v2]}, 
  With[{t = 
     a/Ω}, (Sin[(1 - t) Ω]  v1 + 
      Sin[t  Ω] v2)/Sin[Ω]]]

n = 1000;
v0 = SparseArray[1 -> 1, n]; (* your vector to rotate *)
test = slerpAngle[v0, RandomPoint[Sphere[n]], ArcCos[0.9]];
test . v0 (* 0.9 *)

If you aren't dealing with normalized vectors, then just record the norm, slerp the normalized vector, and then multiply by the original norm at the end.