Generate a random SO(3) matrix.

generate_random_so3_matrix.py

from random import uniform
from math import cos, sin, pi, acos
from sympy import Symbol, I, Matrix, factor, simplify, im, re, expand, true

EPSILON = 0.0001
MAX_ITERATIONS = 1000

a = Symbol('a', real=true)
b = Symbol('b', real=true)
c = Symbol('c', real=true)
d = Symbol('d', real=true)


def newtons_method(y_val):
	f = lambda x : x - 0.5 * sin(2*x) - pi * y_val
	fprime = lambda x : 1 - cos(2 * x) - pi

	# initial value
	x = pi / 2 
	num_iterations = 0

	while( (abs(f(x)) > EPSILON) | (num_iterations > MAX_ITERATIONS) ):
		x += f(x) / fprime(x)
		num_iterations+=1

	return x

def spherical_to_r4_coords(theta1, theta2, theta3):
	x = sin(theta3) * sin(theta2) * cos(theta1)
	y = sin(theta3) * sin(theta2) * sin(theta1)
	z = sin(theta3) * cos(theta2)
	w = cos(theta3)

	return [x, y, z, w] 

def generate_s3_point():
	x, y, z = uniform(0,1), uniform(0,1), uniform(0,1)
	theta1 = 2 * pi * x
	theta2 = acos(1-2*y)
	theta3 = newtons_method(z)

	return theta1, theta2, theta3

def generate_so3_coords():
	x = a + I*b
	z = c + I*d

	# basis of su(2) lie algebra; a real vector space
	e1 = Matrix([[I,0],[0,-I]])
	e2 = Matrix([[0,-1],[1,0]])
	e3 = Matrix([[0,I],[I,0]])

	# generic SU(2) matrix
	g = Matrix([[x, -z.conjugate()],[z, x.conjugate()]])

	# adjoint map from su(2) into su(2)
	ad = lambda X : factor(simplify(expand(g * X * g.H)))

	# represent su(2) matrix as a vector in the basis e1, e2, e3
	SU2Coeff = lambda A : [im(A[0,0]), re(A[1,0]), im(A[1,0])]

	# adjoint representation as a matrix
	B = Matrix([SU2Coeff(ad(e1)), SU2Coeff(ad(e2)), SU2Coeff(ad(e3))]).transpose()
	return B

def generate_random_so3_mat():
	B = generate_so3_coords()
	t1, t2, t3 = generate_s3_point()
	x,y,z,w = spherical_to_r4_coords(t1, t2, t3)
	return B.subs(a, x).subs(b, y).subs(c, z).subs(d, w)