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.

  • $\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$ Commented Sep 26, 2017 at 20:10

11 Answers 11


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.

  • 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.
    Commented 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$ Commented Dec 29, 2020 at 21:42
  • $\begingroup$ This method doesn't seem to generate uniform rotations, as visualized in my answer. Even after applying @AlecJacobson's half-space correction, the result still contains some non-uniformity. My answer implements a few other methods that do seem to work, though: randn_orthobasis, randn_axis ($\psi$ corrected), nbubis ($\psi$ corrected). $\endgroup$ Commented Jan 9 at 9:31
  • $\begingroup$ What I ended up using is here github.com/alecjacobson/gptoolbox/blob/master/matrix/… and doesn't seem to have the same biases as @MateenUlhaq's python code. $\endgroup$ Commented Jan 10 at 14:32
  • $\begingroup$ Ah, with random column sign-flipping (and "optional" determinant correction) demonstrated by @AlecJacobson, qr_full now generates visually uniform rotations. $\endgroup$ Commented Jan 11 at 2:12

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

  • 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$ Commented 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$ Commented 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$ Commented Jul 13, 2013 at 2:13
  • 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$ Commented Jul 13, 2013 at 5:17
  • $\begingroup$ According to Miles 1965 ( 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$ $\endgroup$
    – sega_sai
    Commented Jun 2, 2020 at 22:23

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.


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

  • 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$ Commented 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$ Commented Jul 13, 2013 at 13:30

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

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

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:


Empirical visualizations of various answers.

$\psi$ corrected
$\psi$ corrected
$\psi \sim \mathcal{U}(0, 2\pi)$
$\psi \sim \mathcal{U}(0, 2\pi)$
$(1, 0, 0)$
$(\frac{1}{9}, -\frac{4}{9}, \frac{8}{9})$

The pink point denotes "initial" point. The blue points are the result of rotating the initial point using 5000 randomly generated 3x3 rotation matrices.

import matplotlib.pyplot as plt
import numpy as np

# https://math.stackexchange.com/questions/442418/random-generation-of-rotation-matrices/1602779#1602779
# http://home.lu.lv/~sd20008/papers/essays/Random%20unitary%20[paper].pdf
# https://github.com/alecjacobson/gptoolbox/blob/master/matrix/rand_rotation.m
def qr_full(num_samples=1):
    z = np.random.randn(num_samples, 3, 3)
    q, r = np.linalg.qr(z)
    sign = 2 * (np.diagonal(r, axis1=-2, axis2=-1) >= 0) - 1
    rot = q
    rot *= sign[..., None, :]
    rot[:, 0, :] *= np.linalg.det(rot)[..., None]
    return rot

# https://math.stackexchange.com/questions/442418/random-generation-of-rotation-matrices/1288873#1288873
# https://math.stackexchange.com/questions/44689/how-to-find-a-random-axis-or-unit-vector-in-3d/44701#44701
def randn_orthobasis(num_samples=1):
    z = np.random.randn(num_samples, 3, 3)
    z = z / np.linalg.norm(z, axis=-1, keepdims=True)
    z[:, 0] = np.cross(z[:, 1], z[:, 2], axis=-1)
    z[:, 0] = z[:, 0] / np.linalg.norm(z[:, 0], axis=-1, keepdims=True)
    z[:, 1] = np.cross(z[:, 2], z[:, 0], axis=-1)
    z[:, 1] = z[:, 1] / np.linalg.norm(z[:, 1], axis=-1, keepdims=True)
    return z

# https://math.stackexchange.com/questions/442418/random-generation-of-rotation-matrices/4394036#4394036
# https://math.stackexchange.com/questions/44689/how-to-find-a-random-axis-or-unit-vector-in-3d/44701#44701
def randn_axis(num_samples=1, corrected=True):
    u = np.random.uniform(0, 1, size=num_samples)
    z = np.random.randn(num_samples, 1, 3)
    z = z / np.linalg.norm(z, axis=-1, keepdims=True)

    if corrected:
        t = np.linspace(0, np.pi, 1024)
        cdf_psi = (t - np.sin(t)) / np.pi
        psi = np.interp(u, cdf_psi, t, left=0, right=np.pi)
        psi = 2 * np.pi * u

    return rot3x3_from_axis_angle(z, psi)

# https://math.stackexchange.com/questions/442418/random-generation-of-rotation-matrices/442423#442423
# https://math.stackexchange.com/questions/44689/how-to-find-a-random-axis-or-unit-vector-in-3d/44691#44691
def nbubis(num_samples=1, corrected=True):
    u1 = np.random.uniform(0, 1, size=num_samples)
    u2 = np.random.uniform(0, 1, size=num_samples)
    u3 = np.random.uniform(0, 1, size=num_samples)

    theta = np.arccos(2 * u1 - 1)
    phi = 2 * np.pi * u2
    axis_vector = [
        np.sin(theta) * np.cos(phi),
        np.sin(theta) * np.sin(phi),
    axis_vector = np.stack(axis_vector, axis=1).reshape(-1, 1, 3)

    if corrected:
        t = np.linspace(0, np.pi, 1024)
        cdf_psi = (t - np.sin(t)) / np.pi
        psi = np.interp(u3, cdf_psi, t, left=0, right=np.pi)
        psi = 2 * np.pi * u3

    return rot3x3_from_axis_angle(axis_vector, psi)

# https://math.stackexchange.com/questions/442418/random-generation-of-rotation-matrices/1602779#1602779
def qr_half(num_samples=1):
    z = np.random.randn(num_samples, 3, 3)
    q, r = np.linalg.qr(z)
    return q

# https://en.wikipedia.org/wiki/Rodrigues%27_rotation_formula#Matrix_notation
def rot3x3_from_axis_angle(axis_vector, angle):
    angle = np.atleast_1d(angle)[..., None, None]
    K = np.cross(np.eye(3), axis_vector)
    return np.eye(3) + np.sin(angle) * K + (1 - np.cos(angle)) * (K @ K)

def plot_scatter(pointses, filename, kwargses):
    fig = plt.figure()
    ax = fig.add_subplot(projection="3d", computed_zorder=False)
    for points, kwargs in zip(pointses, kwargses):
        ax.scatter(*np.asarray(points).T, marker=".", **kwargs)
    ax.view_init(elev=45, azim=-45, roll=0)
    ax.set(xlim=(-1, 1), ylim=(-1, 1), zlim=(-1, 1))
    ax.set_aspect("equal", adjustable="box")
    fig.savefig(filename, dpi=300, bbox_inches="tight", pad_inches=0)
    # plt.show()
    # exit()

    "randn_orthobasis": randn_orthobasis,
    "randn_axis": randn_axis,
    "randn_axis_incorrect": lambda **kwargs: randn_axis(corrected=False, **kwargs),
    "nbubis": nbubis,
    "nbubis_incorrect": lambda **kwargs: nbubis(corrected=False, **kwargs),
    "qr_half": qr_half,
    "qr_full": qr_full,

# x is the starting point; y contains various sample rotated points.
# x = np.array([1.0, 0.0, 0.0])
x = np.array([1 / 9, -4 / 9, 8 / 9])
x /= np.linalg.norm(x)  # Normalize to unit vector, just in case.

for name, func in METHODS.items():
    rot = func(num_samples=5000 // (2 if "_half" in name else 1))
    y = rot @ x
        [y, [x]],
        [{"s": 1, "alpha": 0.5}, {"s": 64, "color": "#ff77cc"}],


  • qr_half does not cover the full space of rotations.
  • qr_full corrects qr_half by randomly flipping the signs of the columns of the resulting rotation matrix, while imposing the required constraint $\det T = 1$ so that $T \in \mathcal{SO}(3)$.
  • randn_orthobasis generates two random vectors, and then generates an orthonormal basis from this by taking a bunch of cross products. A more generic orthonormalization method may be used instead for higher dimension matrices (e.g., 4x4, 5x5, ...).
  • randn_axis and nbubis both generate a random 3x3 rotation matrix using a spherically-uniform random axis vector, then a random angle of rotation about this vector. However, choosing $\psi \sim \mathcal{U}(0, 2\pi)$ has a counterintuitive consequence: the initial points are more likely to remain near their starting positions. To repair this, we must make it less likely to choose small angles by choosing $\psi \sim \frac{d}{du} \left( \frac{1}{\pi} (u - \sin u) \right)$, where $u \sim \mathcal{U}(0, \pi)$. In practice, we can use "inverse transform sampling" to sample from the corresponding CDF by simply interpolating over its inverse.
  • $\begingroup$ I haven't looked at this question and its answers in detail in a while. So I have missed your answer here. Thank you very much for this. This was obviously alot of work and to my mind, is very useful for the apparently many people who need to make sense of the various answers here. Thanks alot!!! $\endgroup$ Commented Jul 10 at 1:27

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
  • $\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$ Commented 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$ Commented Jul 13, 2013 at 6:08

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;

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;

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

  • 4
    $\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
    Commented 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$ Commented 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$ Commented Jul 16, 2013 at 22:57

If you have MATLAB, I suggest using the code from What Is a Random Orthogonal Matrix? by the late Professor Nicholas Higham (n in code represents the dimension).

[Q,R] = qr(randn(n));
Q = Q * diag(sign(diag(R)));

If you want to avoid matrix multiplication, you can do it more sufficiently.

[Q,R] = qr(randn(n));
Q(:, diag(R) < 0) = -Q(:, diag(R) < 0);

You must log in to answer this question.

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