40
$\begingroup$

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.

$\endgroup$
2
  • 1
    $\begingroup$ Bob - did you solve this with any code? Would you be open to sharing your code? $\endgroup$
    – jason m
    Sep 26, 2017 at 17:42
  • $\begingroup$ 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. $\endgroup$ Sep 26, 2017 at 20:10

9 Answers 9

20
$\begingroup$

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.

$\endgroup$
3
  • 1
    $\begingroup$ A much better answer for practical purposes than the one marked by OP. $\endgroup$
    – Ben
    Aug 30, 2020 at 11:16
  • 1
    $\begingroup$ 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; $\endgroup$
    – Jan M.
    Sep 22, 2020 at 12:42
  • $\begingroup$ 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); $\endgroup$ Dec 29, 2020 at 21:42
17
$\begingroup$

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]$.

$\endgroup$
9
  • 4
    $\begingroup$ 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. $\endgroup$ Jul 13, 2013 at 2:03
  • 9
    $\begingroup$ @bob.sacamento if you choose the $\theta$ and $\phi$ like I wrote (using $\arccos$), the density will not be higher at the poles. $\endgroup$
    – nbubis
    Jul 13, 2013 at 2:11
  • 6
    $\begingroup$ @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. $\endgroup$
    – nbubis
    Jul 13, 2013 at 2:13
  • 2
    $\begingroup$ yes, I understand your method now. It will give a uniform distribution. Thanks. $\endgroup$ Jul 13, 2013 at 3:46
  • 10
    $\begingroup$ 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+$. $\endgroup$ Jul 13, 2013 at 5:17
15
$\begingroup$

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.

$\endgroup$
1
11
$\begingroup$

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$...)

$\endgroup$
2
  • 2
    $\begingroup$ +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. $\endgroup$ Jul 13, 2013 at 5:21
  • 1
    $\begingroup$ @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. $\endgroup$ Jul 13, 2013 at 13:30
8
$\begingroup$

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)$.

$\endgroup$
3
  • $\begingroup$ I think this has the advantage of working for higher dimensions and not just for 3 dimensions $\endgroup$
    – gota
    Oct 9, 2017 at 20:34
  • 1
    $\begingroup$ Does it? It uses the cross product... $\endgroup$
    – N. McA.
    May 8, 2018 at 14:18
  • 1
    $\begingroup$ 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. $\endgroup$
    – bcmpinc
    May 8, 2018 at 17:21
4
$\begingroup$

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:

$\endgroup$
1
$\begingroup$

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
$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$ Jul 13, 2013 at 2:00
  • $\begingroup$ 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) $\endgroup$ Jul 13, 2013 at 6:08
1
$\begingroup$

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;
}
$\endgroup$
0
$\begingroup$

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).

$\endgroup$
3
  • 3
    $\begingroup$ No. You cannot just use Euler angles generated from a uniform distribution. See Miles 1965 and the link in @FixedPoint's answer. $\endgroup$
    – horchler
    Jul 16, 2013 at 22:41
  • 1
    $\begingroup$ @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. $\endgroup$ Jul 16, 2013 at 22:52
  • $\begingroup$ @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. $\endgroup$ Jul 16, 2013 at 22:57

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .