Random rotation of dimension n

rand_rot.py

import numpy as np
import numpy.linalg as la

def rand_unit(n):
    r = np.random.randn(n)
    return r/la.norm(r)

def householder(unit_v):
    return np.eye(len(unit_v))  - 2*np.outer(unit_v,unit_v)

def random_rotation_matrix(n):
    arr = np.eye(n)
    for dim in range(1,n):
        ref = householder(rand_unit(dim+1))
        arr[:dim+1,:dim+1] = ref@arr[:dim+1,:dim+1]
    #Positive determinant
    if not n%2:
        arr[0] *= -1
    return arr