34

For a current project, I need to generate several $3\times 3$ rotation matrices for input into an algorithm. I thought I might go about this by randomly generating the number of elements needed to define a rotation matrix and then calculating the remaining elements from them. Does anyone know of an algorithm for calculating the remaining elements once a defining set of elements is given? Or does anyone know of a better way to go about this? Thanks.

Ben Grossmann
  • 203,051
  • 12
  • 142
  • 283
bob.sacamento
  • 3,768
  • 2
  • 15
  • 23
  • 1
    Bob - did you solve this with any code? Would you be open to sharing your code? – jason m Sep 26 '17 at 17:42
  • This was from a while back. Might take a while to dig the code out. I pretty much implemented nbubis's answer, except that I took into account the more complicated distribution for $\psi$ that you'll see in the fifth comment on his answer. – bob.sacamento Sep 26 '17 at 20:10

9 Answers9

18

I had this question and I solved it via QR decomposition of 3x3 matrices with normally distributed random numbers. If you have matlab, use: [Q,R] = qr(randn(3));

Q is a uniformly random rotation matrix. For more info, refer to: M. Ozols, “How to generate a random unitary matrix,” 2009.

user2844647
  • 301
  • 2
  • 2
  • A much better answer for practical purposes than the one marked by OP. – Ben Aug 30 '20 at 11:16
  • 1
    Is it possible that this only generates rotations on the half-sphere? Consider (Matlab): `for i=1:1000; [Q,~] = qr(randn(3)); dirs(:,i) = Q*[1 0 0]'; end; scatter3(dirs(1,:), dirs(2,:), dirs(3,:)); axis equal vis3d;` – Jan M. Sep 22 '20 at 12:42
  • for n=3, this produces det(Q)=1, but in general it appears to produce det(Q) = ±1 depending on whether n is odd/even. To fix this and the half-space issue that Jan M. mentions, you could use: `[Q,~] = qr(randn(n));Q(:,1)=Q(:,1)*(2*(rand>0.5)-1);Q(:,2)=det(Q)*Q(:,2);` – Alec Jacobson Dec 29 '20 at 21:42
16

Rotation matrices can be uniquely defined by a vector and a rotation angle. To generate the vector, you can use grandom spherical coordinates $\phi$ and $\theta$. Thus you should first generate random angles by using:

$$\theta = \arccos(2u_1 - 1)$$ $$\phi = 2\pi u_2$$ Where $u_1,u_2$ are uniformly distributed in $[0,1]$. This will give you a vector around which to rotate. Then, randomly decide the amount of rotation $\psi\in[0,2\pi]$.

nbubis
  • 31,733
  • 7
  • 76
  • 133
  • 3
    I like this, with one caveat. I the "density" of $\theta$ and $\phi$ is much higher at the poles than at the equator. Maybe a more truly random distribution would be obtained by randomly selecting $x,y,z$, then computing the corresponding $\theta$ and $\phi$ each time? What do you think? Thanks. – bob.sacamento Jul 13 '13 at 02:03
  • 9
    @bob.sacamento if you choose the $\theta$ and $\phi$ like I wrote (using $\arccos$), the density will **not** be higher at the poles. – nbubis Jul 13 '13 at 02:11
  • 6
    @bob.sacamento, also note that randomly selecting x,y,z, will give you a higher density at the 8 "corner" points, vs. the method I suggested that will give you a uniform distribution on the sphere. – nbubis Jul 13 '13 at 02:13
  • 2
    yes, I understand your method now. It will give a uniform distribution. Thanks. – bob.sacamento Jul 13 '13 at 03:46
  • 10
    Hmm. If uniform means, as I think it should, uniform w.r.t. the Haar measure on SO(3), then according to the link given by Fixed Point the cdf of $\psi$ should be $(\psi-\sin\psi)/\pi, 0\le\psi\le\pi$. IOW you cannot generate $\psi$ uniformly. I didn't check this though. But it does make sense in a way. All the rotations with $\psi\approx0$ are close to each other, and by differentiating you see that the density of $\psi$ goes to zero as $\psi\to0+$. – Jyrki Lahtonen Jul 13 '13 at 05:17
  • According to Miles 1965 ( https://www.jstor.org/stable/2333716) The rotation angle has to be chosen with probability distribution of $\frac{2}{\pi} sin^2(\frac{\phi}{2})$ where $0<\phi<\pi$ – sega_sai Jun 02 '20 at 22:23
  • @sega_sai The interval of $\phi$ should be $0\leq\phi<2\pi$, to make the integral being 1. – chaosink Aug 06 '20 at 08:22
  • @chaosink you are wrong. $\int_0^\infty sin^2(x/2) dx= \pi/2$ – sega_sai Aug 09 '20 at 04:43
  • @sega_sai $\displaystyle\int_0^{2\pi}\cfrac 2 \pi sin^2(\cfrac \psi 2)d\psi=\cfrac{\psi-sin\psi}{2\pi}\Bigg|_0^{2\pi}=1$ – chaosink Aug 11 '20 at 02:16
13

I have had the same exact problem myself a while ago so I point you to this which says it very succinctly with plenty of references.

Fixed Point
  • 7,524
  • 3
  • 28
  • 46
  • the wikipedia references are covered pretty cleanly in this paper https://statweb.stanford.edu/~cgates/PERSI/papers/subgroup-rand-var.pdf – MeowBlingBling May 05 '21 at 01:06
11

This is an interesting question! If the "random matrix" is being used for any sort of Monte-Carlo testing, it is important that a sequence of "random matrices" produced by whatever (as von Neumann noted, putting us in a state of sin if we claim too much about it) pseudo-random device, be at least demonstrably equi-distributed in the rotation group $SO(3)$.

I have no serious operational quibble with other answers, which are certainly thoughtful and productive, but/and I might object that they are ad-hoc, so offering no a-priori promise of any genuine pseudo-randomness.

If it does matter to have a more-certifiable pseudo-randomness, the following device lends itself more to proof, for random 3-D rotations. Use the fact that Hamiltonian $\mathbb H$ quaternions give 3D rotations in at least one way, namely, identify $\mathbb R^3$ with purely imaginary quaternions, and let $g\in \mathbb H^\times$ act on purely-imaginary quaternions $x$ by $g\cdot x=gxg^{-1}$.

In that set-up, it's not too hard to prove various "pseudo-randomness" (equi-distribution) properties.

(Edit: For example, choose elements of $\mathbb H$ pseudo-randomly in some ball of radius $R$ according to some volume-equidistribution "rule", and let $R\rightarrow +\infty$...)

paul garrett
  • 46,394
  • 4
  • 79
  • 149
  • 2
    +1 for getting uniform w.r.t. the Haar measure. The unit quaternion suffice (they are a double cover of SO(3)), so no need to let $R\to\infty$ :-) You do need to generate points in 4D with a spherically symmetric distribution. – Jyrki Lahtonen Jul 13 '13 at 05:21
  • 1
    @JyrkiLahtonen, one still might want to choose integer quaternions uniformly distributed in a ball of size $R$, and divide by $R$ to make unit-norm ones... to effect the choice of random rotations in a discretized way. – paul garrett Jul 13 '13 at 13:30
8

A rotation matrix R is the same as an orthonormal basis that preserves orientation ($\det(R)=1$). Hence, to create a uniform distributed random rotation matrix, we need to pick three orthogonal random unit vectors, make sure that the orientation is correct and concatenate them into a matrix.

First pick two random unit vectors $u$ and $v$. For numerical stability, ensure that $|u\cdot v| < 0.99$. Now project $v$ onto the plane perpendicular to $u$, i.e. subtract $(v\cdot u)u$ from $v$ and normalize the result. The last axis is fixed by $u$ and $v$, hence we can compute $w$ using the cross product: $w = u \times v$. The resulting rotation matrix is $R=(u\ v\ w)$.

bcmpinc
  • 219
  • 2
  • 3
  • I think this has the advantage of working for higher dimensions and not just for 3 dimensions – gota Oct 09 '17 at 20:34
  • 1
    Does it? It uses the cross product... – N. McA. May 08 '18 at 14:18
  • 1
    The cross product is not strictly necessary. You can pick $w$ at random and project it onto the space perpendicular to the plane $u,v$ and normalize the resulting vector. This can be repeated for more vectors, though at some point numerical stability will become an issue. If the resulting matrix has $\det(R) = -1$, invert the last vector. Also, picking unit vectors in higher dimensional spaces is more difficult. – bcmpinc May 08 '18 at 17:21
4

There are many ways to do this. One method is described as follows,

To generate uniformly distributed random rotations of a unit sphere, first perform a random rotation about the vertical axis, then rotate the north pole to a random position.

from Fast Random Rotation Matrices by James Arvo

Code realization:

Syrtis Major
  • 183
  • 6
1

We can take advantage of the Rodrigues' rotation formula to generate rotation matrices.

Image of the Rodrigues' rotation formula

Image of the bracket notation

With this, you only need to generate a random 3x1 unit vector (as the rotation axis) and specify a rotation angle. This formula will transform them into a valid rotation matrix.

MATLAB example:

function R = rot(w, theta)
  bw = [0, -w(3), w(2); w(3), 0, -w(1); -w(2), w(1), 0];
  R = eye(3) + sin(theta)*bw + (1-cos(theta))*bw*bw;
end

w = rand(3,1)
w = w/norm(w)
R = rot(w, 3.14)

C++ example:

// w: the unit vector indicating the rotation axis
// theta: the rotation angle in radian
Eigen::Matrix3d MatrixExp3 (Eigen::Vector3d w, float theta){
  Eigen::Matrix3d bw, R;
  bw << 0, -w(2), w(1), w(2), 0, -w(0), -w(1), w(0), 0;
  R << Eigen::Matrix3d::Identity() + std::sin(theta)*bw + (1-std::cos(theta))*bw*bw;
  return R;
}

int main() {
  std::srand((unsigned int) time(0));
  Eigen::Vector3d w = Eigen::Vector3d::Random();
  Eigen::Matrix3d R = MatrixExp3(w.normalized(), 3.14f);
  std::cout << R << std::endl;
}
0

Here's one way: use yaw, pitch, and roll. Generate three angles at random: $\alpha,\beta,\gamma$. Generate your random rotation as

$$ R = R_x(\gamma) \, R_y(\beta) \, R_z(\alpha)\ $$

Where $R_z(\alpha)$ is the rotation matrix about the $z$-axis, given by $$ R_z = \begin{bmatrix} \cos\alpha & -\sin\alpha & 0\\ \sin\alpha & \cos\alpha & 0\\ 0 & 0 & 1 \end{bmatrix} $$

and the other two are similarly defined (for their respective axes and angles).

Ben Grossmann
  • 203,051
  • 12
  • 142
  • 283
  • 3
    No. You cannot just use Euler angles generated from a uniform distribution. See [Miles 1965](http://dx.doi.org/10.2307%2F2333716) and the link in @FixedPoint's answer. – horchler Jul 16 '13 at 22:41
  • 1
    @horchler I thought that this method was sufficient given the parameters of the question: this is certainly a way of randomly generating a rotation matrix, and a consicely explained way at that. I never claimed that all rotations were of equal likelihood; I only said that plugging in the Euler angles will give you a rotation matrix. The advantage of this method is that all rotations can be generated, if perhaps with varying likelihoods. – Ben Grossmann Jul 16 '13 at 22:52
  • @horchler Also strictly speaking, these are not the Euler angles, though I'm sure they are isomorphic in some regard. Thank you for the reference though, it's interesting that higher dimensional rotations have this counterintuitive property. – Ben Grossmann Jul 16 '13 at 22:57
0

What I often do is, to generate a random matrix and column-rotate it to some shape, for instance triangular shape. keeping track of the rotations gives then a random-rotation-matrix.
In my matrix-program MatMate this read simply:

 u = randomu(3,3)   // generate a 3x3 matrix with uniform distributed entries
 t = gettrans ( u, "tri") // get the rotation-matrix, which rotates u to "tri"angular shape

 // check
 chk = u * t   // matrix chk is triangular
Gottfried Helms
  • 32,738
  • 3
  • 60
  • 134
  • Very interesting. I am not familiar with the technique. Could you point me to a proof that the resulting matrix will always be a rotation matrix. Not that I doubt your word. I'm just wanting to edify myself. Thanks. – bob.sacamento Jul 13 '13 at 02:00
  • Well, such a proof is not really needed - it is just, that from the random general matrix the angles are determined, which are taken in the answer of @omnom themselves as random. Only the ***distribution*** of the angles are different: while in Omnom's answer they are taken from approximate equidistributed (uniform) points from the surface of a sphere, my procedure assumes points uniformly equidistributed in the full volume of a cube. Therefor the distribution of the random rotation matrices is different (and I cannot know whch method would serve your purposes better) – Gottfried Helms Jul 13 '13 at 06:08