Trying out code from "Using Eigendecomposition to Convert Rotations and Interpolate Operations". Found some numerical stability issues. https://algassert.com/quantum/2016/01/10/eigendecomposition-for-rotation-and-interpolation.html

interpolation_stability.py

import numpy as np
from scipy.stats import unitary_group

def eigenterpolate(U0, U1, s):
	"""Interpolates between two matrices."""
	return U0 * eigenpow(U0.H * U1, s)

def eigenpow(M, t):
	"""Raises a matrix to a power."""
	return eigenlift(lambda b: b**t, M)

def eigenlift(f, M):
	"""Lifts a numeric function to apply it to a matrix."""
	w, v = np.linalg.eig(M)
	T = np.mat(np.zeros(M.shape, np.complex128))
	for i in range(len(w)):
		eigen_val = w[i]
		eigen_vec = np.mat(v[:, i])
		eigen_mat = np.dot(eigen_vec, eigen_vec.H)
		T += f(eigen_val) * eigen_mat
	return T

def is_unitary(m):
	return np.allclose(np.eye(len(m)), m.dot(m.T.conj()))
	#return np.allclose(np.eye(m.shape[0]), np.asmatrix(m).H * m)

def rotation_change(u, rate):
	return eigenterpolate(np.asmatrix(u), unitary_group.rvs(2), rate)

# Generate a random unitary matrix
u = unitary_group.rvs(2)
# Iterate for a hundred generations
for i in range(100):
	# Try to slightly change the rotation
	u_ = rotation_change(u, 0.001)
	# If unitarity is broken, retry changing the rotation
	for j in range(1000):
		if is_unitary(u_): break
		print('Retrying interpolation...')
		u_ = rotation_change(u, 0.001)
	# Update the unitary matrix if we found a small change that still keeps it
	# unitary, or exit if we failed to do so
	u = u_ if is_unitary(u_) else exit()
	# Print information on the matrix
	print(u)
	print(f'Is unitary: {is_unitary(u)}')