An important component of the Cartan KAK decomposition for 2 qubit operations is to diagonalize a 4x4 unitary matrix using orthogonal (not unitary, purely real orthogonal) matrices. That is to say, given unitary U find orthogonal A and B such that
A*U*B is diagonal. (Actually, the orthogonal matrices are supposed to be special orthogonal but that's easily fixed.)
Writing code to do this correctly (nevermind quickly) is a giant pain. Is there a method included with numpy that could be used to do most of the heavy lifting?
For example, this problem reduces to simultaneously diagonalizing
imag(U). As a wild guess I tested if the svd of
real(U) + random.random()*imag(U) would give a result that happened to work. Numpy does give orthogonal matrices in this situation, but they don't always diagonalize the original U unfortunately.