# Diagonalize a unitary matrix with orthogonal matrices using numpy

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 real(U) and 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.

• Can't you just use the eigenvalue decomposition? Commented Sep 11, 2019 at 2:57
• @nicoguaro I think the issue here is (complex) unitary vs. (real) orthogonal. The eigendecomposition of a complex matrix returns complex $A$ and $B$, in general. Commented Sep 11, 2019 at 8:05

Have you tried the QZ decomposition on real(U) and imag(U)? In general it returns AA and BB upper triangular rather than diagonal, but I wonder if a stroke of luck happens here (exactly like the Schur decomposition can be proved to always return a diagonal $$T$$ when run on a normal matrix, rather than an upper triangular one).

(EDIT: yes, on second thought this "stroke of luck" should happen, because AA + i*BB is both unitary and upper triangular, so it must be diagonal.)

(second EDIT: a quick test on Matlab (which should use the exact same LAPACK libraries for this decomposition) confirms the stroke of luck.)

import numpy as np
import scipy.linalg
def diagonalize_unitary_using_two_orthogonals(u):
"""Decomposes u into L @ np.diag(D) @ R.T where L and R are real orthogonal.
"""
diag_r, diag_i, left, right = scipy.linalg.qz(np.real(u), np.imag(u))
diag = np.diagonal(diag_r) + np.diagonal(diag_i) * 1j
return left, diag, right

• I checked this for myself and it does appear to work. Awesome! I don't suppose there's a nice way to also guarantee the results are special orthogonal instead of just orthogonal? It's not too hard to fix after the fact by checking the determinants and inverting a column or row, though. Commented Sep 11, 2019 at 20:35
• @CraigGidney In the Lapack code, those orthogonal matrices are assembled as products of Householder reflectors (det=-1) and Givens rotations (det=1), so it's an easy modification to keep track of their determinants along the code and swap a sign at the end. However you need to modify the Lapack code to do that; I don't think there is a way to do it directly from the Python wrapper without dirtying your hands with Fortran 77 code. Commented Sep 11, 2019 at 21:00
• This fails on some matrices that are almost equal to the identity. For example here is a case that fails to diagonalize: np.array([[(1.0000000000000002+0j), 0j, 0j, 4.266421588589642e-17j], [0j, (1.0000000000000002+0j), (-0-4.266421588589642e-17j), 0j], [0j, 4.266421588589642e-17j, (1.0000000000000002+0j), 0j], [(-0-4.266421588589642e-17j), 0j, 0j, (1.0000000000000002+0j)]]) Commented Sep 11, 2019 at 22:47
• @CraigGidney Uh, I see the problem. Actually, to keep the output real, the real version of QZ is allowed to return matrices that have 2x2 diagonal blocks, instead of being upper triangular (or diagonal in our case). That's how they can handle real matrices with complex eigenvalues, for instance. I'll be thinking if there is an easy workaround --- all we have left to do now is handling 2x2 subproblems. But certainly that's less elegant than I was imagining. Commented Sep 12, 2019 at 6:25
• The only error cases I found were cases that resulted in a permutation matrix instead of a diagonal matrix. Are you aware of a way to get a non-permutation matrix as the "diagonal" result? Commented Sep 12, 2019 at 20:40