I want to sample uniformly at random from all $n\times n$ orthogonal matrices $M$, subject to two requirements:
$M$ is an orthogonal matrix. Loosely speaking, it represents a rotation (or rotation with reflection).
$Mx$ is at angle $\le \theta$ from $x$, for all $x$, where $\theta$ is a fixed constant. Equivalently, $x^\top Mx \le c x^\top x$ for all $x$, where $c=\cos^{-1} \theta$ is a fixed constant.
Loosely speaking, these matrices represent rotation matrices that rotate by at most a bounded amount.
Is there an algorithm or method to sample uniformly at random from the set of such matrices?
I'm working in high dimensions, where $n$ is large (say 1000; not 2 or 3). My application is probably too esoteric to be of general interest; it relates to generating new unit vectors that are "similar" to existing unit vectors, with similarity measured by angle / cosine distance. Rejection sampling seems ineffective in high-dimensional spaces, as the proportion of orthogonal matrices that satisfy the second constraint is exponentially small, so it would require exponentially many draws from $O(n)$ to find a single such matrix.