4
$\begingroup$

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.

Improve this question
$\endgroup$
2
  • $\begingroup$ Can't you just use the eigenvalue decomposition? $\endgroup$
    – nicoguaro ♦
    Sep 11, 2019 at 2:57
  • 1
    $\begingroup$ @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. $\endgroup$
    – Federico Poloni
    Sep 11, 2019 at 8:05

1 Answer 1

Reset to default
6
$\begingroup$

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
Improve this answer
$\endgroup$
6
  • $\begingroup$ 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. $\endgroup$
    – Craig Gidney
    Sep 11, 2019 at 20:35
  • $\begingroup$ @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. $\endgroup$
    – Federico Poloni
    Sep 11, 2019 at 21:00
  • $\begingroup$ 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)]]) $\endgroup$
    – Craig Gidney
    Sep 11, 2019 at 22:47
  • 1
    $\begingroup$ @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. $\endgroup$
    – Federico Poloni
    Sep 12, 2019 at 6:25
  • $\begingroup$ 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? $\endgroup$
    – Craig Gidney
    Sep 12, 2019 at 20:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.