Probabilistically correct SO(3) sampling

[brentyi]

Motivation and Context

The current uniform sampling procedure for SO(3) works by:

1. Uniformly sampling in log/tangent space.
2. Applying the exponential map.

This would produce a uniform output if the exponential operation were linear, but it's unfortunately not.

As intuition for why there might be something funky afoot:

• If we consider a subset of the domain of exp_map(omega) where l2_norm(omega) = 2 * pi, all inputs maps to a singleton.
• If we consider a subset of the domain of exp_map(omega) where l2_norm(omega) = pi /2, the function is injective, and produces a unique rotation for each unique tangent vector.

How Has This Been Tested

Unit tests still pass!

The implementation is also battle-tested:

• This is a near-identical PR as this 2019 one to Sophus: Probabilistically correct uniform sampling for SO(3) strasdat/Sophus#233
• It's a port from the jaxlie implementation.

Types of changes

• Docs change / refactoring / dependency upgrade
• Bug fix (non-breaking change which fixes an issue)
• New feature (non-breaking change which adds functionality)
• Breaking change (fix or feature that would cause existing functionality to change)

Checklist

• My code follows the code style of this project.
• My change requires a change to the documentation.
• I have updated the documentation accordingly.
• I have read the CONTRIBUTING document.
• I have completed my CLA (see CONTRIBUTING)
• I have added tests to cover my changes.
• All new and existing tests passed.

[fantaosha]

LGTM! Thanks for this nice feature!

[luisenp]

Thanks a lot for the nice contribution @brentyi! I have a question for you and @fantaosha. Do we need similar procedures for our other sampling methods? In thinking in particular of randn() and also those in SO2 class.

[fantaosha]

Thanks a lot for the nice contribution @brentyi! I have a question for you and @fantaosha. Do we need similar procedures for our other sampling methods? In thinking in particular of randn() and also those in SO2 class.

No, we don't. This problem only happens to SO3 and SE3. Since SE3.rand() depends on SO3.rand(), we only need to refactor SO3.rand()

[brentyi]

Appreciate the fast review!!

Can confirm that SO2.rand also looks fine, this is much simpler than sampling on SO(3).

I'm wondering for SO2.randn and SO3.randn: these look correct if the goal is to sample from some zero-mean unit-covariance tangent-space Gaussians, but what is that actually used for? With Gaussians in Euclidean spaces randn ~ N(0, I) is a reasonable abstraction because a mean/covariance can be imposed by scaling shifting/scaling sampled vectors, but this doesn't seem to extend in a principled way to the SO2/SO3 cases.

[mhmukadam]

Thanks for the contribution, @brentyi.

[mhmukadam]

for SO2.randn and SO3.randn ... what is that actually used for?

Good question. For now I think we mostly use it in unit tests but I imagine it could be used to sample say a random (group) trajectory as input to some optimization problem. You are right the current Gaussian sampling is in tangent space. Maybe @fantaosha can expand on any other general use case where what we have currently is insufficient.

[fantaosha]

for SO2.randn and SO3.randn ... what is that actually used for?

Good question. For now I think we mostly use it in unit tests but I imagine it could be used to sample say a random (group) trajectory as input to some optimization problem. You are right the current Gaussian sampling is in tangent space. Maybe @fantaosha can expand on any other general use case where what we have currently is insufficient.

Hi @brentyi, as @mhmukadam said, rand() and randn() originally referred to the random generators used to sample the tangent space of the identity. Honestly, we didn't consider any specific geometric meaning of uniform or normal distribution on Lie groups. However, we are happy to refactor randn() as what you have contributed to rand(). Thanks!

[empet]

Searching the web to find out where the Shoemake method of sampling uniformly from the group of unit quaternions is implemented, I stumpled upon this thread, "Probabilistically correct SO(3) sampling". On math.stackexchange there is an unanswered question related to the correctness of the Shoemake algorithm that was implemented here by @brentyi. We'd be grateful if you could explain what led you to admit that the algorithm samples uniformly a unit quaternion (and as a consequence an element in SO(3)).

[brentyi]

Hi!

Two links that might be of interest to you are:

• Planning Algorithms (LaValle 2006): https://web.archive.org/web/20211105205926/http://planning.cs.uiuc.edu/node198.html
• UnitRandom() in Eigen: https://eigen.tuxfamily.org/dox/Quaternion_8h_source.html

[empet]

The presentation from the web archive link refers to the Shoemake paper, [883], which is incriminated on math.stackexchange.
In addition, it is explained in a blog style what Hopf fibration is, without shedding more light on Shoemake's approach.

EIGEN just implements operations on quaternions, including the Shoemake method of uniform sampling from SO(3).
Unfortunately a rigurous explanation of this method is still lacking. Just read the quoted text from math.stackexchange and try to understand it and perform the product of the two quaternions involved there. Comparing with the Shoemake result you'll see that it was considered that the sum of two independent uniform random variables is also uniformly distributed, which is FALSE.
Thank you for trying to give links you considered to be more convincing.