def rvs(dim=3):
    random_state = np.random
H = np.eye(dim)
D = np.ones((dim,))
for n in range(1, dim):
    x = random_state.normal(size=(dim - n + 1,))
D[n - 1] = np.sign(x[0])
x[0] -= D[n - 1] * np.sqrt((x * x).sum())
Hx = np.eye(dim - n + 1) - 2.0 * np.outer(x, x) / (x * x).sum()
mat = np.eye(dim)
mat[n - 1:, n - 1:] = Hx
H = np.dot(H, mat)
D[-1] = (-1) ** (1 - dim % 2) * D.prod()
H = (D * H.T).T
return H
How to create random orthonormal matrix in python numpy