> 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. James Arvo, Fast Random Rotation Matrices

rand_rotmat.py

def rand_rotmat(n):
    # http://www.realtimerendering.com/resources/GraphicsGems/gemsiii/rand_rotation.c
    import numpy as np

    theta, phi, z = np.random.rand(3, n) * 2
    theta *= np.pi
    phi *= np.pi

    r = np.sqrt(z)
    Vx = np.sin(phi) * r
    Vy = np.cos(phi) * r
    Vz = np.sqrt(2.0 - z)

    st = np.sin(theta)
    ct = np.cos(theta)
    Sx = Vx * ct - Vy * st
    Sy = Vx * st + Vy * ct

    M = np.empty((n, 3, 3), 'float')

    M[:, 0, 0] = Vx * Sx - ct
    M[:, 0, 1] = Vx * Sy - st
    M[:, 0, 2] = Vx * Vz

    M[:, 1, 0] = Vy * Sx + st
    M[:, 1, 1] = Vy * Sy - ct
    M[:, 1, 2] = Vy * Vz

    M[:, 2, 0] = Vz * Sx
    M[:, 2, 1] = Vz * Sy
    M[:, 2, 2] = 1.0 - z
    return M

test.py

from pylab import *

a = rand_rotmat(10000).T
print (a[0]*a[1]).sum(0).max()

t = arctan2(a[0, 0], a[0, 1])
hist(t, 41, normed=True);
axhline(0.5/pi)

t = arctan2(a[0, 0], (a[0, 1]**2+a[0, 2]**2)**0.5)
hist(t, 41, normed=True);
x = linspace(-pi, pi, 41)/2
plot(x, cos(x)/2)


a, b = rand(2, 10, 3)
m = rand_rotmat(10)
a2 = matmul(m, a[..., newaxis]).squeeze(-1)
b2 = matmul(m, b[..., newaxis]).squeeze(-1)
assert allclose((a*b).sum(-1), (a2*b2).sum(-1))