MathOverflow is a question and answer site for professional mathematicians. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

An answer to another question derived a formula for the volume of a delta-ball in $O(n)$. I am wondering if there is a (constructive) way to draw samples uniformly at random from such a region.

For my purposes, an acceptance-rejection scheme is not viable. I am trying to study the behavior of an integral over $O(n)$ of a function that is very small outside of certain regions, so accepting samples that lie within the ball is no different (from a computational perspective) from sampling from $O(n)$ in its entirety.

Alternatively, any other non-uniform sampling scheme on $O(n)$ might be of interest, especially if it concentrates more samples around diagonal matrices $\text{diag}(\pm1,\dots,\pm1)$.

share|cite|improve this question
For small $\delta$, sampling uniformly from a $\delta$-ball in the (linear) space of antisymmetric matrices and exponentiating ought to give an approximately uniform sample. – Mark Meckes Dec 19 '11 at 23:45
If you tell what your function is, somebody may come up with a better integration scheme than Monte-Carlo, not excluding the possibility of the exact value or good asymptotics. – fedja Dec 20 '11 at 0:21
@Mark Meckes: I guess you could even do an importance sampling based on the Jacobian of the exponential map. It doesn't seem so hard to compute the Jacobian. Another interesting note (though not of particular relevance): Persi and Saloff-Coste mentioned in their paper on the subgroup algorithm that there is a bijection on a set of measure 1 from SO(n) to open sets of the Euclidean space, which looks very similar to the stereographical projection. It's invented by H. Weyl. – John Jiang Dec 20 '11 at 3:10
up vote 5 down vote accepted

In case anyone else comes across this: the model I ended up using is by León, Massé, and Rivest in the Journal of Multivariate Analysis (see here). They give a distribution on the space of skew-symmetric matrices that gives an arbitrarily concentrated distribution on the orthogonal group: that is, for dispersion parameter $\kappa=0$ the resulting distribution is uniform, and as $\kappa\rightarrow\infty$ the distribution approaches a constant distribution.

EDIT: The limit is a Gaussian distribution (analogous to the degrees of freedom of a t distribution approaching $\infty$), not constant.

Thanks again to those who commented/answered as these suggestions put me on the path towards the method I ultimately found.

share|cite|improve this answer

I have not proven the following to work, but I've done similar calculations when debugging software in my 2009-10 work with Matt Hastings. If these seem fast enough, it should not be too much work to validate them.

(1) I would take small, random perturbations of the identity matrix and then do four or five iterations of Newtons method to get the polar part of the resulting matrix. I offer this suggestion under the assumption that you can tolerate some samples falling outside the delta ball.

By Newton's method I mean replacing $X$ by $\frac{1}{2}\left( X + (X^{-1})^\mathrm{T}\right) $ where $T$ is the transpose. This is a simplification of what is suggested by N. J. Higham for computing the polar part of an invertible matrix.

The number of iterations needed depends on delta. I have no idea how big your matrices are. If that are bigger than 200 by 200 you might need to use a package that has a well parallelized matrix inversion routine. If you are looking at smaller matrices, you might just go ahead and apply 10 iterations.

(2) If you cannot abide by matrices outside the delta ball, then generate a random diagonal orthogonal $D$ that is close to $I$ and a random orthogonal $W$ and multiply to $U = WDW^\mathrm{T}$ as the desired random orthogonal, now known to be in the delta ball. I am assuming the operator norm, which is costly to compute by the default in pure math, I think.

You need to generate the random unitary by taking polar decomposition of a random matrix. Now you need to read up on Newton's method for for computing the polar part, as you may deal with badly conditioned matrices.

I am not sure what distribution you want for the diagonal orthogonal matrix. So I am hoping (1) will work for you.

share|cite|improve this answer
Thanks very much for the response. I am currently experimenting with implementing these ideas and if it works out then I'll accept (I honestly don't have any idea if it will or not, but I'm hopeful). – bnaul Dec 28 '11 at 20:48

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.