seeking to simplify `find-diagonalizer-in-e-basis`

[stylewarning]

EDIT This issue was previously about using generalized Schur to re-implement F-D-I-E-B, but it turned out to not be e a working approach. As such, I'm changing this issue into a discussion about how we might approach doing so.

I'm mostly interpreting some things @kilimanjaro told me (so credit goes to him), though errors I may have below are my own.

find-diagonalizer-in-e-basis aims to diagonalize a symmetric unitary matrix gammag = u u^T in terms of an orthogonal matrix of eigenvectors. Normally, all that the spectral theorem gives you (in this context) is that you can do this in terms of a unitary matrix of eigenvectors. In this case, however, the real and imaginary parts of gammag commute, so it's sufficient to simultaneously diagonalize them. They are real, symmetric matrices, their eigenvectors will be real and will give an orthogonal matrix

To simultaneously diagonalize a pair of commuting matrices, it's not quite enough to just compute eigenvectors and eigenvalues of one and then hope that this works for the other (consider that the identity matrix, which commutes with anything, but we need to write its eigenspaces as spanned by eigenvectors of some other matrix.) The current QUILC approach is to try to diagonalize a linear combination of the real and imaginary parts. Since there could be some relationships between them that we don't know, the current code tries to pick a random combination, and just repeats until we get something that works

This problem has a standard solution that linear algebra libraries implement, called the generalized Schur decomposition. LAPACK documents it here

A Python implementation (which uses NumPy's qz) supplied by @kilimanjaro is here:

def orthogonal_decomposition(U):
    """
    Given unitary U, decompose UU^T into O @ np.diag(d) @ O.T where O is special orthogonal.
    """
    # Generalized Schur decomposition writes
    #    Re[U] = L D_r R^T,  Im[U] = L D_i R^T
    # In general, D_r, D_i are upper triangular, but for unitary U they end up diagonal
    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
    if np.linalg.det(left) < 0:
        left[:,0] *= -1
    return left, diag*diag

• [x] If we want to implement this in QUILC, we'll need to get qz into MAGICL.

[braised-babbage]

I was wrong about this. The real QZ decomposition A,B = Q*AA*Z^h, Q*BB*Z^h given by dgges (cf. here) can and does result in 2x2 blocks in AA on occasion. In fact, this appears to be quite sensitive to small perturbations of A,B. For example, consider the following

(defvar u1 (magicl:from-list '(#C(0.7071067690849304d0 0.0d0) #C(-0.7071067690849304d0 0.0d0) #C(0.0d0 0.0d0)
 #C(0.0d0 0.0d0) #C(0.7071067690849304d0 0.0d0) #C(0.7071067690849304d0 0.0d0)
 #C(0.0d0 0.0d0) #C(0.0d0 0.0d0) #C(0.0d0 0.0d0) #C(0.0d0 0.0d0)
 #C(0.7071067690849304d0 0.0d0) #C(-0.7071067690849304d0 0.0d0) #C(0.0d0 0.0d0)
 #C(0.0d0 0.0d0) #C(0.7071067690849304d0 0.0d0) #C(0.7071067690849304d0 0.0d0)) '(4 4)))
U1
CL-QUIL> (defvar u2 (magicl:from-list '(#C(0.7071067811865476d0 -5.551115123125783d-17)
 #C(-0.7071067811865477d0 1.1102230246251565d-16) #C(0.0d0 0.0d0)
 #C(0.0d0 0.0d0) #C(0.7071067811865477d0 -1.1102230246251565d-16)
 #C(0.7071067811865476d0 -5.551115123125783d-17) #C(0.0d0 0.0d0)
 #C(0.0d0 0.0d0) #C(0.0d0 0.0d0) #C(0.0d0 0.0d0)
 #C(0.7071067811865476d0 5.551115123125783d-17)
 #C(-0.7071067811865477d0 -1.1102230246251565d-16) #C(0.0d0 0.0d0)
 #C(0.0d0 0.0d0) #C(0.7071067811865477d0 1.1102230246251565d-16)
 #C(0.7071067811865476d0 5.551115123125783d-17)) '(4 4)))
U2
CL-QUIL> (magicl:norm (magicl:reshape (magicl:.- u1 u2) '(16)))
3.4228542365326185d-8
CL-QUIL> (defun qzri (mat) (magicl:qz (magicl:.realpart mat) (magicl:.imagpart mat)))
QZRI
CL-QUIL> (qzri u1)
#<MAGICL:MATRIX/DOUBLE-FLOAT (4x4):
   1.000     0.000     0.000     0.000
   0.000     1.000     0.000     0.000
   0.000     0.000     1.000     0.000
   0.000     0.000     0.000     1.000>
#<MAGICL:MATRIX/DOUBLE-FLOAT (4x4):
   0.000     0.000     0.000     0.000
   0.000     0.000     0.000     0.000
   0.000     0.000     0.000     0.000
   0.000     0.000     0.000     0.000>
#<MAGICL:MATRIX/DOUBLE-FLOAT (4x4):
   1.000     0.000     0.000     0.000
  -0.000     1.000     0.000     0.000
  -0.000     0.000     1.000     0.000
  -0.000     0.000     0.000     1.000>
#<MAGICL:MATRIX/DOUBLE-FLOAT (4x4):
   0.707     0.707     0.000     0.000
  -0.707     0.707     0.000     0.000
   0.000     0.000     0.707     0.707
   0.000     0.000    -0.707     0.707>
CL-QUIL> (qzri u2)
#<MAGICL:MATRIX/DOUBLE-FLOAT (4x4):
  -0.949     0.316     0.000     0.000
  -0.316    -0.949    -0.000    -0.000
   0.000     0.000     0.949    -0.316
   0.000     0.000     0.316     0.949>
#<MAGICL:MATRIX/DOUBLE-FLOAT (4x4):
   0.000    -0.000     0.000     0.000
   0.000     0.000    -0.000    -0.000
   0.000     0.000     0.000    -0.000
   0.000     0.000     0.000     0.000>
#<MAGICL:MATRIX/DOUBLE-FLOAT (4x4):
  -0.230    -0.973     0.000     0.000
  -0.973     0.230     0.000     0.000
   0.000     0.000    -0.230    -0.973
   0.000     0.000    -0.973     0.230>
#<MAGICL:MATRIX/DOUBLE-FLOAT (4x4):
   0.973     0.230    -0.000    -0.000
   0.230    -0.973    -0.000    -0.000
   0.000    -0.000    -0.973    -0.230
   0.000    -0.000    -0.230     0.973>

Note the significant divergence between the factorizations, despite the small difference between u1 and u2. In this case, u2 gets 2x2 blocks in the "diagonal" part of the factorization.

You can see how deceptive the above can be. For some reason I was under the impression that, because UU^T is symmetric, that when you compute Q*AA*Z^h and Q*BB*Z^h, then UU^T = Q*(AA+i*BB)^2*Q^T would have the (AA+i*BB)^2 turn out to be diagonal due to symmetry. This is also not true in general.

The Python code I had (by the way, qz was not even my idea, I found it from some stackexchange post somewhere) was half baked and not thoroughly tested. It did work on the examples that I tried, but that's not enough in this case.

The broader approach of this orthogonal factorization is that we are trying to take advantage of the commutation of the real and imaginary parts of a unitary matrix U, by instead simultaneously diagonalizing these. According to random strangers this is a hard problem, but there may be some algorithms (e.g. the mentioned one using "Jacobi angles") which can tackle this. On the other hand, there seem to be some strings attached, and these algorithms are not standard lapack routines.

[braised-babbage]

For what it's worth, I spent a couple of hours scouring LAPACK to see if I could spot an easy routines for what we want, but to no success. It's still not clear to me whether there's something clever that we can do here that I'm missing.

[genos]

qutip has a simdiag function, though others have noted some numerical instability.

[genos]

There's also NUMERICAL METHODS FOR SIMULTANEOUS DIAGONALIZATION, which this rando Matlab file claims to implement.

[genos]

Since U U^T is symmetric, Golub and Van Loan (pg. 500 in the 4th edition) may have what we need (though I don't know if B = Im{U U^T} would be positive definite):

Copy-pasta of the above screenshot as text

Algorithm 8.7.1 Given A = A^T ∈ IR^{n×n} and B = B^T ∈ IR{n×n} with B positive definite, the following algorithm computes a nonsingular X such that X^T A X = diag(a1 , . . . , an ) and X^T B X = I_n. Compute the Cholesky factorization B = G G^T using Algorithm 4.2.2. Compute C = G^{−1} A G^{−T}. Use the symmetric QR algorithm to compute the Schur decomposition Q^T C Q = diag(a1, . . . , an). Set X = G^{−T} Q.

[braised-babbage]

Since U U^T is symmetric, Golub and Van Loan (pg. 500 in the 4th edition) may have what we need (though I don't know if B = Im{U U^T} would be positive definite):

I think U U^T can be basically any symmetric unitary, so we're not going to generally have positive definiteness of real or imaginary parts. For example, UU^T = np.diag([1, -1, 1j, -1j]]) or any conjugate of this.

[genos]

Sure, but I think we can modify the above approach since we don't need X^T B X to be I. Replacing the Cholesky decomposition with an eigenvalue one (np.linalg.eig uses "the _geev LAPACK routines"), I think we may be in business:

import numpy as np
from scipy.stats import unitary_group
from rich import print
from rich.progress import track
import typer

def decomp_uut(u: np.ndarray) -> np.ndarray:
    uut = u @ u.T
    a, b = uut.real, uut.imag
    _, g = np.linalg.eig(b)
    g_1 = np.linalg.pinv(g)
    c = g_1 @ a @ g_1.T
    _, v = np.linalg.eig(c)
    return g_1.T @ v.T

def assert_is_almost_diag(x: np.ndarray):
    assert (d := np.abs(x - np.diag(np.diag(x))).max()) < 1e-8, f"Max off diag: {d}"

def test_decomp_uut(
    iterations: int = typer.Option(100_000, help="Number of tests to run"),
    dim: int = typer.Option(10, help="Dimension of square matrices"),
    seed: int = typer.Option(8675309, help="PRNG seed"),
):
    print(f"[magenta]Running {iterations:,} test(s)[/magenta]")
    rng = unitary_group(dim=dim, seed=seed)
    for _ in track(range(iterations), description="[cyan]Testing…[/cyan]"):
        u = rng.rvs()
        uut = u @ u.T
        a, b = uut.real, uut.imag
        x = decomp_uut(u)
        assert_is_almost_diag(x.T @ a @ x)
        assert_is_almost_diag(x.T @ b @ x)
    print("[bold][green]Passed![/green][/bold]")

if __name__ == "__main__":
    typer.run(test_decomp_uut)

[genos]

Untested, but I think this should do in Lisp:

(defun decomp-uut (u)
  (let* ((uut (magicl:@ u (magicl:transpose u)))
         (a (magicl:.realpart uut))
         (b (magicl:.imagpart uut)))
    (multiple-value-bind (_ g) (magicl:eig b)
      (declare (ignore _))
      (let* ((g-inv (magicl:inv g))
             (g-inv-transpose (magicl:transpose g-inv))
             (c (magicl:@ g-inv a g-inv-transpose)))
        (multiple-value-bind (_ v) (magicl:eig c)
          (declare (ignore _))
          (magicl:@ g-inv-transpose (magicl:transpose v)))))))

[stylewarning]

@genos Looking into this now.

[stylewarning]

@genos I implemented this in https://github.com/quil-lang/quilc/pull/850; I don't have more specific feedback but I'm finding the following.

First, for random unitaries, it seems to work. I'm essentially running your code verbatim, except I'm orthogonalizing after and ensuring determinant = 1.

CL-QUIL> (loop :repeat 100000 :do
           (let ((m (random-unitary 4)))
             (diagonalizer-in-e-basis m))
           :finally (print 'success))

SUCCESS 

This includes passing a couple math tests. However, when running within QUILC, I get errors, namely:

X^T (UU^T) X not diagonal!

X =
#<MATRIX/COMPLEX-DOUBLE-FLOAT (4x4):
   0.680 + 0.000j     0.056 + 0.000j    -0.288 + 0.000j    -0.672 + 0.000j
   0.438 + 0.000j     0.752 + 0.000j     0.332 + 0.000j     0.363 + 0.000j
   0.228 + 0.000j    -0.428 + 0.000j     0.858 + 0.000j    -0.173 + 0.000j
   0.542 + 0.000j    -0.498 + 0.000j    -0.267 + 0.000j     0.622 + 0.000j>

U =
#<MATRIX/COMPLEX-DOUBLE-FLOAT (4x4):
  -0.354 - 0.354j     0.354 - 0.354j     0.354 - 0.354j    -0.354 - 0.354j
  -0.354 + 0.354j    -0.354 - 0.354j     0.354 + 0.354j     0.354 - 0.354j
  -0.354 + 0.354j     0.354 + 0.354j    -0.354 - 0.354j     0.354 - 0.354j
  -0.354 - 0.354j    -0.354 + 0.354j    -0.354 + 0.354j    -0.354 - 0.354j>

UU^T =
#<MATRIX/COMPLEX-DOUBLE-FLOAT (4x4):
   0.000 + 0.000j     0.000 - 0.000j     0.000 + 0.000j    -0.000 + 1.000j
   0.000 - 0.000j     0.000 + 0.000j     0.000 - 1.000j     0.000 + 0.000j
   0.000 + 0.000j     0.000 - 1.000j     0.000 - 0.000j    -0.000 - 0.000j
  -0.000 + 1.000j     0.000 + 0.000j    -0.000 - 0.000j     0.000 - 0.000j>

X^T(UU^T)X =
#<MATRIX/COMPLEX-DOUBLE-FLOAT (4x4):
   0.000 + 0.538j     0.000 - 0.292j     0.000 - 0.789j    -0.000 + 0.051j
   0.000 - 0.292j    -0.000 + 0.588j     0.000 - 0.374j    -0.000 + 0.655j
   0.000 - 0.789j     0.000 - 0.374j     0.000 - 0.415j    -0.000 - 0.254j
  -0.000 + 0.051j    -0.000 + 0.655j    -0.000 - 0.254j     0.000 - 0.710j>

Original M such that U = E^T M E is
#<MATRIX/COMPLEX-DOUBLE-FLOAT (4x4):
  -0.707 - 0.707j    -0.000 - 0.000j    -0.000 + 0.000j    -0.000 + 0.000j
  -0.000 + 0.000j     0.000 + 0.000j    -0.000 - 0.000j    -0.707 - 0.707j
   0.000 + 0.000j    -0.000 + 0.000j    -0.707 - 0.707j    -0.000 + 0.000j
   0.000 - 0.000j    -0.707 - 0.707j    -0.000 + 0.000j    -0.000 - 0.000j>

These are similar errors to what we were getting with @kilimanjaro's approach. It seems that low-dimensional subsets of the unitary group are particularly troublesome.

[stylewarning]

Continuing the last message, we see that $M$ is essentially the matrix for $\mathtt{CNOT}\;0\;1$ with an extra factor. If we plug this in directly, we get some cleaner results.

CL-QUIL> (magicl:from-list '(1 0 0 0
                          0 0 0 1
                          0 0 1 0
                          0 1 0 0)
                           '(4 4)
                            :type '(complex double-float))
#<MAGICL:MATRIX/COMPLEX-DOUBLE-FLOAT (4x4):
   1.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j
   0.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j     1.000 + 0.000j
   0.000 + 0.000j     0.000 + 0.000j     1.000 + 0.000j     0.000 + 0.000j
   0.000 + 0.000j     1.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j>
CL-QUIL> (find-diagonalizer-in-e-basis *)
X^T (UU^T) X not diagonal!

X =
#<MATRIX/COMPLEX-DOUBLE-FLOAT (4x4):
  -0.707 + 0.000j     0.000 + 0.000j     0.000 + 0.000j     0.707 + 0.000j
  -0.707 + 0.000j     0.000 + 0.000j     0.000 + 0.000j    -0.707 + 0.000j
   0.000 + 0.000j     0.707 + 0.000j    -0.707 + 0.000j     0.000 + 0.000j
   0.000 + 0.000j     0.707 + 0.000j     0.707 + 0.000j     0.000 + 0.000j>

U =
#<MATRIX/COMPLEX-DOUBLE-FLOAT (4x4):
   0.500 + 0.000j     0.000 + 0.500j     0.000 + 0.500j     0.500 + 0.000j
   0.000 - 0.500j     0.500 + 0.000j    -0.500 + 0.000j     0.000 + 0.500j
   0.000 - 0.500j    -0.500 + 0.000j     0.500 + 0.000j     0.000 + 0.500j
   0.500 + 0.000j     0.000 - 0.500j     0.000 - 0.500j     0.500 + 0.000j>

UU^T =
#<MATRIX/COMPLEX-DOUBLE-FLOAT (4x4):
   0.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j     1.000 + 0.000j
   0.000 + 0.000j     0.000 + 0.000j    -1.000 + 0.000j     0.000 + 0.000j
   0.000 + 0.000j    -1.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j
   1.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j>

X^T(UU^T)X =
#<MATRIX/COMPLEX-DOUBLE-FLOAT (4x4):
   0.000 + 0.000j     0.000 + 0.000j    -1.000 + 0.000j     0.000 + 0.000j
   0.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j     1.000 + 0.000j
  -1.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j
   0.000 + 0.000j     1.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j>

Original M such that U = E^T M E is
#<MATRIX/COMPLEX-DOUBLE-FLOAT (4x4):
   1.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j
   0.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j     1.000 + 0.000j
   0.000 + 0.000j     0.000 + 0.000j     1.000 + 0.000j     0.000 + 0.000j
   0.000 + 0.000j     1.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j>
   [Condition of type SIMPLE-ERROR]

[stylewarning]

Disregard, had a typo.

Continuing the last comment, it may be tempting to think that this is because we were simultaneously diagonalizing a real matrix with a zero matrix, which seems like bunk, but we'll get the same problem even if we pass a complex diagonal matrix in, such as $\mathtt{CPHASE}(\pi/2)\;1\;0 = \mathrm{diag}(1,1,1,i)$, which (1) ought to not need diagonalizing, but (2) is in turn is diagonalizing $\mathrm{diag}(1,1,1,0)$ and $\mathrm{diag}(0,0,0,1)$:

CL-QUIL> (magicl:from-list '(1 0 0 0
                          0 1 0 1
                          0 0 1 0
                          0 0 0 #C(0 1))
                           '(4 4)
                            :type '(complex double-float))
#<MAGICL:MATRIX/COMPLEX-DOUBLE-FLOAT (4x4):
   1.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j
   0.000 + 0.000j     1.000 + 0.000j     0.000 + 0.000j     1.000 + 0.000j
   0.000 + 0.000j     0.000 + 0.000j     1.000 + 0.000j     0.000 + 0.000j
   0.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j     0.000 + 1.000j>
CL-QUIL> (find-diagonalizer-in-e-basis *)
X^T (UU^T) X not diagonal!

X =
#<MATRIX/COMPLEX-DOUBLE-FLOAT (4x4):
  -0.946 + 0.000j    -0.041 + 0.000j     0.062 + 0.000j    -0.315 + 0.000j
  -0.297 + 0.000j     0.497 + 0.000j    -0.235 + 0.000j     0.781 + 0.000j
  -0.125 + 0.000j    -0.769 + 0.000j     0.320 + 0.000j     0.539 + 0.000j
   0.031 + 0.000j     0.399 + 0.000j     0.916 + 0.000j     0.033 + 0.000j>

U =
#<MATRIX/COMPLEX-DOUBLE-FLOAT (4x4):
   0.500 + 0.500j     0.500 + 0.500j     0.000 + 0.000j     0.000 + 0.000j
  -0.500 - 0.500j     0.500 + 0.500j     0.000 + 0.000j     0.000 + 0.000j
   0.000 - 0.500j    -0.500 + 0.000j     1.000 + 0.000j     0.000 + 0.000j
   0.500 + 0.000j     0.000 - 0.500j     0.000 + 0.000j     1.000 + 0.000j>

UU^T =
#<MATRIX/COMPLEX-DOUBLE-FLOAT (4x4):
   0.000 + 1.000j     0.000 + 0.000j     0.000 - 0.500j     0.500 + 0.000j
   0.000 + 0.000j     0.000 + 1.000j    -0.500 + 0.000j     0.000 - 0.500j
   0.000 - 0.500j    -0.500 + 0.000j     1.000 + 0.000j     0.000 + 0.000j
   0.500 + 0.000j     0.000 - 0.500j     0.000 + 0.000j     1.000 + 0.000j>

X^T(UU^T)X =
#<MATRIX/COMPLEX-DOUBLE-FLOAT (4x4):
  -0.050 + 0.875j    -0.164 - 0.424j    -0.411 + 0.306j     0.042 + 0.294j
  -0.164 - 0.424j     1.117 + 0.019j    -0.057 - 0.270j    -0.298 + 0.128j
  -0.411 + 0.306j    -0.057 - 0.270j     1.073 + 0.254j    -0.002 - 0.523j
   0.042 + 0.294j    -0.298 + 0.128j    -0.002 - 0.523j    -0.140 + 0.852j>

Original M such that U = E^T M E is
#<MATRIX/COMPLEX-DOUBLE-FLOAT (4x4):
   1.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j
   0.000 + 0.000j     1.000 + 0.000j     0.000 + 0.000j     1.000 + 0.000j
   0.000 + 0.000j     0.000 + 0.000j     1.000 + 0.000j     0.000 + 0.000j
   0.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j     0.000 + 1.000j>
   [Condition of type SIMPLE-ERROR]

Maybe the real and imaginary parts each must be non-singular for this to work?

[genos]

1. Thanks @stylewarning and @kilimanjaro for working so hard on this with me!
2. For cases such as the above where one or both parts are singular, is there a better approach? And how well does MAGICL support checking for singularity/looking at condition numbers, etc.?
3. OMG GitHub markdown supports $\LaTeX$ now!?

[stylewarning]

1. Of course!

2. It's pretty scrappy. We could try to improve things here.

3. Yes! Though it's still a little rough (e.g. sometimes needing to double escape, sometimes not).

[genos]

Continuing with the $\mathtt{CPHASE}(\pi/2) 1 0$ example, I think there's a typo in your definition @stylewarning; the second row has an extra 1 in it. Switching the definition to be $\mathtt{diag}(1, 1, 1, i)$ looks to work in both Python and Lisp:

from icecream import ic
import numpy as np

ic(cphase := np.diag([1, 1, 1, 1j]))  # cphase(π / 2) 1 0
a, b = cphase.real, cphase.imag
ic(a, b)
b_vals, g = np.linalg.eig(b)
ic(b_vals, g)
ic(g_inv := np.linalg.pinv(g))
ic(c := g_inv @ a @ g_inv.T)
c_vals, v = np.linalg.eig(c)
ic(c_vals, v)
ic(x := g_inv.T @ v.T)
ic(x.T @ a @ x)
ic(x.T @ b @ x)

python cphase.py
ic| cphase := np.diag([1, 1, 1, 1j]): array([[1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],
                                             [0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j],
                                             [0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j],
                                             [0.+0.j, 0.+0.j, 0.+0.j, 0.+1.j]])
ic| a: array([[1., 0., 0., 0.],
              [0., 1., 0., 0.],
              [0., 0., 1., 0.],
              [0., 0., 0., 0.]])
    b: array([[0., 0., 0., 0.],
              [0., 0., 0., 0.],
              [0., 0., 0., 0.],
              [0., 0., 0., 1.]])
ic| b_vals: array([0., 0., 0., 1.])
    g: array([[1., 0., 0., 0.],
              [0., 1., 0., 0.],
              [0., 0., 1., 0.],
              [0., 0., 0., 1.]])
ic| g_inv := np.linalg.pinv(g): array([[1., 0., 0., 0.],
                                       [0., 1., 0., 0.],
                                       [0., 0., 1., 0.],
                                       [0., 0., 0., 1.]])
ic| c := g_inv @ a @ g_inv.T: array([[1., 0., 0., 0.],
                                     [0., 1., 0., 0.],
                                     [0., 0., 1., 0.],
                                     [0., 0., 0., 0.]])
ic| c_vals: array([1., 1., 1., 0.])
    v: array([[1., 0., 0., 0.],
              [0., 1., 0., 0.],
              [0., 0., 1., 0.],
              [0., 0., 0., 1.]])
ic| x := g_inv.T @ v.T: array([[1., 0., 0., 0.],
                               [0., 1., 0., 0.],
                               [0., 0., 1., 0.],
                               [0., 0., 0., 1.]])
ic| x.T @ a @ x: array([[1., 0., 0., 0.],
                        [0., 1., 0., 0.],
                        [0., 0., 1., 0.],
                        [0., 0., 0., 0.]])
ic| x.T @ b @ x: array([[0., 0., 0., 0.],
                        [0., 0., 0., 0.],
                        [0., 0., 0., 0.],
                        [0., 0., 0., 1.]])

(ql:quickload :magicl)

(let* ((cphase (magicl:from-list '(1 0 0 0
                                   0 1 0 0
                                   0 0 1 0
                                   0 0 0 #C(0 1))
                                 '(4 4)
                                 :type '(complex double-float)))
       (a (magicl:.realpart cphase))
       (b (magicl:.imagpart cphase)))
  (multiple-value-bind (b-vals g) (magicl:eig b)
    (let* ((g-inv (magicl:inv g))
           (g-inv-transpose (magicl:transpose g-inv))
           (c (magicl:@ g-inv a g-inv-transpose)))
      (multiple-value-bind (c-vals v) (magicl:eig c)
        (let ((x (magicl:@ g-inv-transpose (magicl:transpose v))))
          (loop :for (name value)
                :in (list `(a ,a)
                          `(b ,b)
                          `(b-vals ,b-vals)
                          `(g ,g)
                          `(g-inv ,g-inv)
                          `(c ,c)
                          `(c-vals ,c-vals)
                          `(v ,v)
                          `(x ,x)
                          `(x^T-a-x ,(magicl:@ (magicl:transpose x) a x))
                          `(x^T-b-x ,(magicl:@ (magicl:transpose x) b x)))
                :do (format t "~A: ~A~%" name value)))))))

sbcl --non-interactive --load cphase.lisp
This is SBCL 2.2.6, an implementation of ANSI Common Lisp.
More information about SBCL is available at <http://www.sbcl.org/>.

SBCL is free software, provided as is, with absolutely no warranty.
It is mostly in the public domain; some portions are provided under
BSD-style licenses.  See the CREDITS and COPYING files in the
distribution for more information.
To load "magicl":
  Load 1 ASDF system:
    magicl
; Loading "magicl"
...............
A: #<MATRIX/DOUBLE-FLOAT (4x4):
      1.000     0.000     0.000     0.000
      0.000     1.000     0.000     0.000
      0.000     0.000     1.000     0.000
      0.000     0.000     0.000     0.000>
B: #<MATRIX/DOUBLE-FLOAT (4x4):
      0.000     0.000     0.000     0.000
      0.000     0.000     0.000     0.000
      0.000     0.000     0.000     0.000
      0.000     0.000     0.000     1.000>
B-VALS: (0.0d0 0.0d0 0.0d0 1.0d0)
G: #<MATRIX/DOUBLE-FLOAT (4x4):
      1.000     0.000     0.000     0.000
      0.000     1.000     0.000     0.000
      0.000     0.000     1.000     0.000
      0.000     0.000     0.000     1.000>
G-INV: #<MATRIX/DOUBLE-FLOAT (4x4):
          1.000     0.000     0.000    -0.000
          0.000     1.000     0.000    -0.000
          0.000     0.000     1.000    -0.000
          0.000     0.000     0.000     1.000>
C: #<MATRIX/DOUBLE-FLOAT (4x4):
      1.000     0.000     0.000     0.000
      0.000     1.000     0.000     0.000
      0.000     0.000     1.000     0.000
      0.000     0.000     0.000     0.000>
C-VALS: (1.0d0 1.0d0 1.0d0 0.0d0)
V: #<MATRIX/DOUBLE-FLOAT (4x4):
      1.000     0.000     0.000     0.000
      0.000     1.000     0.000     0.000
      0.000     0.000     1.000     0.000
      0.000     0.000     0.000     1.000>
X: #<MATRIX/DOUBLE-FLOAT (4x4):
      1.000     0.000     0.000     0.000
      0.000     1.000     0.000     0.000
      0.000     0.000     1.000     0.000
      0.000     0.000     0.000     1.000>
X^T-A-X: #<MATRIX/DOUBLE-FLOAT (4x4):
            1.000     0.000     0.000     0.000
            0.000     1.000     0.000     0.000
            0.000     0.000     1.000     0.000
            0.000     0.000     0.000     0.000>
X^T-B-X: #<MATRIX/DOUBLE-FLOAT (4x4):
            0.000     0.000     0.000     0.000
            0.000     0.000     0.000     0.000
            0.000     0.000     0.000     0.000
            0.000     0.000     0.000     1.000>

[genos]

The $\mathtt{CNOT}$ example still fails, as taking b_vals, g = np.linalg.eig(np.zeros((2, 2))) ensures that g is the identity matrix, so $g^{-1}ag^{-T}$ is still $a$ i.e. $\mathtt{CNOT}$, which isn't diagonal. Maybe we need only check that the input matrix isn't in fact real?

[genos]

Apologies for getting Python all over your QUILC, but I'm still quicker there.

from icecream import ic
import numpy as np

def decompose(u: np.ndarray) -> np.ndarray:
    if np.isreal(u).all():
        _, x = np.linalg.eig(u)
    else:
        a, b = u.real, u.imag
        _, g = np.linalg.eig(b)
        g_inv = np.linalg.pinv(g)
        c = g_inv @ a @ g_inv.T
        _, v = np.linalg.eig(c)
        x = g_inv.T @ v.T
    return x

def is_almost_diag(x: np.ndarray) -> bool:
    return np.abs(x - np.diag(np.diag(x))).max() < 1e-8

cnot = np.array([[1, 0, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0], [0, 1, 0, 0]])
cphase = np.diag([1, 1, 1, 1j])
for u in [cnot, cphase]:
    ic(u)
    a, b = u.real, u.imag
    ic(x := decompose(u))
    ic(x.T @ a @ x)
    ic(is_almost_diag(x.T @ a @ x))
    ic(x.T @ b @ x)
    ic(is_almost_diag(x.T @ b @ x))

python real_or_complex.py
ic| u: array([[1, 0, 0, 0],
              [0, 0, 0, 1],
              [0, 0, 1, 0],
              [0, 1, 0, 0]])
ic| x := decompose(u): array([[ 0.        ,  0.        ,  1.        ,  0.        ],
                              [ 0.70710678,  0.70710678,  0.        ,  0.        ],
                              [ 0.        ,  0.        ,  0.        ,  1.        ],
                              [ 0.70710678, -0.70710678,  0.        ,  0.        ]])
ic| x.T @ a @ x: array([[ 1.,  0.,  0.,  0.],
                        [ 0., -1.,  0.,  0.],
                        [ 0.,  0.,  1.,  0.],
                        [ 0.,  0.,  0.,  1.]])
ic| is_almost_diag(x.T @ a @ x): True
ic| x.T @ b @ x: array([[0., 0., 0., 0.],
                        [0., 0., 0., 0.],
                        [0., 0., 0., 0.],
                        [0., 0., 0., 0.]])
ic| is_almost_diag(x.T @ b @ x): True
ic| u: array([[1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],
              [0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j],
              [0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j],
              [0.+0.j, 0.+0.j, 0.+0.j, 0.+1.j]])
ic| x := decompose(u): array([[1., 0., 0., 0.],
                              [0., 1., 0., 0.],
                              [0., 0., 1., 0.],
                              [0., 0., 0., 1.]])
ic| x.T @ a @ x: array([[1., 0., 0., 0.],
                        [0., 1., 0., 0.],
                        [0., 0., 1., 0.],
                        [0., 0., 0., 0.]])
ic| is_almost_diag(x.T @ a @ x): True
ic| x.T @ b @ x: array([[0., 0., 0., 0.],
                        [0., 0., 0., 0.],
                        [0., 0., 0., 0.],
                        [0., 0., 0., 1.]])
ic| is_almost_diag(x.T @ b @ x): True

[genos]

Hacky Lisp attempt

(ql:quickload :magicl)

(defun real->complex (m)
  "Convert a real matrix M to a complex one."
  (let ((cm (magicl:zeros
              (magicl:shape m)
              :type `(complex ,(magicl:element-type m)))))
    (magicl::map-to #'complex m cm)
    cm))

(defun decompose (u)
  (let* ((a (magicl:.realpart u))
         (b (magicl:.imagpart u))
         (abs-max-b (reduce #'max (magicl::storage (magicl:map #'abs b))))
         (x (if (< abs-max-b 1e-8)
              (multiple-value-bind (_ x) (magicl:eig u)
                (declare (ignore _))
                x)
              (multiple-value-bind (_ g) (magicl:eig b)
                (declare (ignore _))
                (let* ((g-inv (real->complex (magicl:inv g)))
                       (g-inv-transpose (magicl:transpose g-inv))
                       (c (magicl:@ g-inv (real->complex a) g-inv-transpose)))
                  (multiple-value-bind (_ v) (magicl:eig c)
                    (declare (ignore _))
                    (magicl:@ g-inv-transpose (magicl:transpose v))))))))
    x))

(let ((cnot (magicl:from-list '(1 0 0 0
                                0 0 0 1
                                0 0 1 0
                                0 1 0 0)
                              '(4 4)
                              :type '(complex double-float)))
      (cphase (magicl:from-list '(1 0 0 0
                                  0 1 0 0
                                  0 0 1 0
                                  0 0 0 #C(0 1))
                                '(4 4)
                                :type '(complex double-float))))
  (loop :for (name unitary)
        :in (list `(cnot ,cnot)
                  `(cphase ,cphase))
        :do (progn
              (format t "~A: ~A~%" name unitary)
              (let* ((a (magicl:.realpart unitary))
                     (b (magicl:.imagpart unitary))
                     (x (decompose unitary))
                     (x^t (magicl:transpose x))
                     (x^t-a-x (magicl:@ x^t (real->complex a) x))
                     (x^t-b-x (magicl:@ x^t (real->complex b) x)))
                (loop :for (k v)
                      :in (list `(a ,a)
                                `(b ,b)
                                `(x ,x)
                                `(x^t ,x^t)
                                `(x^t-a-x ,x^t-a-x)
                                `(x^t-b-x ,x^t-b-x))
                      :do (format t "~A: ~A~%" k v))))))

sbcl --non-interactive --load real_or_complex.lisp
This is SBCL 2.2.6, an implementation of ANSI Common Lisp.
More information about SBCL is available at <http://www.sbcl.org/>.

SBCL is free software, provided as is, with absolutely no warranty.
It is mostly in the public domain; some portions are provided under
BSD-style licenses.  See the CREDITS and COPYING files in the
distribution for more information.
To load "magicl":
  Load 1 ASDF system:
    magicl
; Loading "magicl"
...............
CNOT: #<MATRIX/COMPLEX-DOUBLE-FLOAT (4x4):
         1.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j
         0.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j     1.000 + 0.000j
         0.000 + 0.000j     0.000 + 0.000j     1.000 + 0.000j     0.000 + 0.000j
         0.000 + 0.000j     1.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j>
A: #<MATRIX/DOUBLE-FLOAT (4x4):
      1.000     0.000     0.000     0.000
      0.000     0.000     0.000     1.000
      0.000     0.000     1.000     0.000
      0.000     1.000     0.000     0.000>
B: #<MATRIX/DOUBLE-FLOAT (4x4):
      0.000     0.000     0.000     0.000
      0.000     0.000     0.000     0.000
      0.000     0.000     0.000     0.000
      0.000     0.000     0.000     0.000>
X: #<MATRIX/COMPLEX-DOUBLE-FLOAT (4x4):
     -0.000 + 0.000j     0.000 + 0.000j     1.000 + 0.000j     0.000 + 0.000j
      0.707 + 0.000j     0.707 + 0.000j     0.000 + 0.000j     0.000 + 0.000j
     -0.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j     1.000 + 0.000j
      0.707 + 0.000j    -0.707 + 0.000j     0.000 + 0.000j     0.000 + 0.000j>
X^T: #<MATRIX/COMPLEX-DOUBLE-FLOAT (4x4):
       -0.000 + 0.000j     0.707 + 0.000j    -0.000 + 0.000j     0.707 + 0.000j
        0.000 + 0.000j     0.707 + 0.000j     0.000 + 0.000j    -0.707 + 0.000j
        1.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j
        0.000 + 0.000j     0.000 + 0.000j     1.000 + 0.000j     0.000 + 0.000j>
X^T-A-X: #<MATRIX/COMPLEX-DOUBLE-FLOAT (4x4):
            1.000 + 0.000j    -0.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j
           -0.000 + 0.000j    -1.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j
            0.000 + 0.000j     0.000 + 0.000j     1.000 + 0.000j     0.000 + 0.000j
            0.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j     1.000 + 0.000j>
X^T-B-X: #<MATRIX/COMPLEX-DOUBLE-FLOAT (4x4):
            0.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j
            0.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j
            0.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j
            0.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j>
CPHASE: #<MATRIX/COMPLEX-DOUBLE-FLOAT (4x4):
           1.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j
           0.000 + 0.000j     1.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j
           0.000 + 0.000j     0.000 + 0.000j     1.000 + 0.000j     0.000 + 0.000j
           0.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j     0.000 + 1.000j>
A: #<MATRIX/DOUBLE-FLOAT (4x4):
      1.000     0.000     0.000     0.000
      0.000     1.000     0.000     0.000
      0.000     0.000     1.000     0.000
      0.000     0.000     0.000     0.000>
B: #<MATRIX/DOUBLE-FLOAT (4x4):
      0.000     0.000     0.000     0.000
      0.000     0.000     0.000     0.000
      0.000     0.000     0.000     0.000
      0.000     0.000     0.000     1.000>
X: #<MATRIX/COMPLEX-DOUBLE-FLOAT (4x4):
      1.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j
      0.000 + 0.000j     1.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j
      0.000 + 0.000j     0.000 + 0.000j     1.000 + 0.000j     0.000 + 0.000j
      0.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j     1.000 + 0.000j>
X^T: #<MATRIX/COMPLEX-DOUBLE-FLOAT (4x4):
        1.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j
        0.000 + 0.000j     1.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j
        0.000 + 0.000j     0.000 + 0.000j     1.000 + 0.000j     0.000 + 0.000j
        0.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j     1.000 + 0.000j>
X^T-A-X: #<MATRIX/COMPLEX-DOUBLE-FLOAT (4x4):
            1.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j
            0.000 + 0.000j     1.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j
            0.000 + 0.000j     0.000 + 0.000j     1.000 + 0.000j     0.000 + 0.000j
            0.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j>
X^T-B-X: #<MATRIX/COMPLEX-DOUBLE-FLOAT (4x4):
            0.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j
            0.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j
            0.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j
            0.000 + 0.000j     0.000 + 0.000j     0.000 + 0.000j     1.000 + 0.000j>

[stylewarning]

Python is OK, and good catch on the CPHASE typo. My bad. I'll look into your proposed change.

[genos]

I suppose we’d probably need a second special case for if the matrix is purely imaginary.

[stylewarning]

I handled the real=0 and imag=0 cases, and this is the next failure I get:

X^T (UU^T) X not diagonal!

X =
#<MATRIX/COMPLEX-DOUBLE-FLOAT (4x4):
   0.259 + 0.000j     0.966 + 0.000j     0.000 + 0.000j     0.000 + 0.000j
   0.000 + 0.000j     0.000 + 0.000j    -0.259 + 0.000j     0.966 + 0.000j
  -0.966 + 0.000j     0.259 + 0.000j     0.000 + 0.000j     0.000 + 0.000j
   0.000 + 0.000j     0.000 + 0.000j     0.966 + 0.000j     0.259 + 0.000j>

U =
#<MATRIX/COMPLEX-DOUBLE-FLOAT (4x4):
   0.924 + 0.000j     0.000 + 0.000j     0.000 - 0.383j     0.000 + 0.000j
   0.000 + 0.000j     0.924 + 0.000j     0.000 + 0.000j     0.000 - 0.383j
   0.000 - 0.383j     0.000 + 0.000j     0.924 + 0.000j     0.000 + 0.000j
   0.000 + 0.000j     0.000 - 0.383j     0.000 + 0.000j     0.924 + 0.000j>

UU^T =
#<MATRIX/COMPLEX-DOUBLE-FLOAT (4x4):
   0.707 + 0.000j     0.000 + 0.000j     0.000 - 0.707j     0.000 + 0.000j
   0.000 + 0.000j     0.707 + 0.000j     0.000 + 0.000j     0.000 - 0.707j
   0.000 - 0.707j     0.000 + 0.000j     0.707 + 0.000j     0.000 + 0.000j
   0.000 + 0.000j     0.000 - 0.707j     0.000 + 0.000j     0.707 + 0.000j>

X^T(UU^T)X =
#<MATRIX/COMPLEX-DOUBLE-FLOAT (4x4):
   0.707 + 0.354j    -0.000 + 0.612j     0.000 + 0.000j     0.000 + 0.000j
  -0.000 + 0.612j     0.707 - 0.354j     0.000 + 0.000j     0.000 + 0.000j
   0.000 + 0.000j     0.000 + 0.000j     0.707 + 0.354j     0.000 - 0.612j
   0.000 + 0.000j     0.000 + 0.000j     0.000 - 0.612j     0.707 - 0.354j>

Original M such that U = E^T M E is
#<MATRIX/COMPLEX-DOUBLE-FLOAT (4x4):
   0.924 + 0.000j     0.000 + 0.000j    -0.383 + 0.000j     0.000 + 0.000j
   0.000 + 0.000j     0.924 + 0.000j     0.000 + 0.000j     0.383 + 0.000j
   0.383 + 0.000j     0.000 + 0.000j     0.924 + 0.000j     0.000 + 0.000j
   0.000 + 0.000j    -0.383 + 0.000j     0.000 + 0.000j     0.924 + 0.000j>
   [Condition of type SIMPLE-ERROR]

[stylewarning]

In this case $UU^T = I\otimes\frac{1}{\sqrt{2}}\begin{pmatrix}1 & -i \\ -i & 1\end{pmatrix}$.

[genos]

Sorry for so much whack-a-mole on this one 😞

I think the issue here is that $a = Im\{UU^T\}$ is already diagonal, in which case we can just take the eigendecomp of $b$.

from icecream import ic
import numpy as np

ic(uut := np.kron(np.eye(2), 1 / np.sqrt(2) * np.array([[1, -1j], [-1j, 1]])))
a, b = uut.real, uut.imag
ic(a, b)
b_vals, g = np.linalg.eig(b)
ic(b_vals, g)
ic(x := np.linalg.pinv(g))
ic(x.T @ a @ x)
ic(x.T @ b @ x)

python decomp.py
ic| uut := np.kron(np.eye(2), 1 / np.sqrt(2) * np.array([[1, -1j], [-1j, 1]])): array([[0.70710678+0.j        , 0.        -0.70710678j,
                                                                                        0.        +0.j        , 0.        +0.j        ],
                                                                                       [0.        -0.70710678j, 0.70710678+0.j        ,
                                                                                        0.        +0.j        , 0.        +0.j        ],
                                                                                       [0.        +0.j        , 0.        +0.j        ,
                                                                                        0.70710678+0.j        , 0.        -0.70710678j],
                                                                                       [0.        +0.j        , 0.        +0.j        ,
                                                                                        0.        -0.70710678j, 0.70710678+0.j        ]])
ic| a: array([[0.70710678, 0.        , 0.        , 0.        ],
              [0.        , 0.70710678, 0.        , 0.        ],
              [0.        , 0.        , 0.70710678, 0.        ],
              [0.        , 0.        , 0.        , 0.70710678]])
    b: array([[ 0.        , -0.70710678,  0.        ,  0.        ],
              [-0.70710678,  0.        ,  0.        ,  0.        ],
              [ 0.        ,  0.        ,  0.        , -0.70710678],
              [ 0.        ,  0.        , -0.70710678,  0.        ]])
ic| b_vals: array([ 0.70710678, -0.70710678,  0.70710678, -0.70710678])
    g: array([[ 0.70710678,  0.70710678,  0.        ,  0.        ],
              [-0.70710678,  0.70710678,  0.        ,  0.        ],
              [ 0.        ,  0.        ,  0.70710678,  0.70710678],
              [ 0.        ,  0.        , -0.70710678,  0.70710678]])
ic| x := np.linalg.pinv(g): array([[ 0.70710678, -0.70710678,  0.        ,  0.        ],
                                   [ 0.70710678,  0.70710678,  0.        ,  0.        ],
                                   [ 0.        ,  0.        ,  0.70710678, -0.70710678],
                                   [ 0.        ,  0.        ,  0.70710678,  0.70710678]])
ic| x.T @ a @ x: array([[7.07106781e-01, 1.66533454e-16, 0.00000000e+00, 0.00000000e+00],
                        [1.89526925e-16, 7.07106781e-01, 0.00000000e+00, 0.00000000e+00],
                        [0.00000000e+00, 0.00000000e+00, 7.07106781e-01, 1.66533454e-16],
                        [0.00000000e+00, 0.00000000e+00, 1.89526925e-16, 7.07106781e-01]])
ic| x.T @ b @ x: array([[-7.07106781e-01, -3.25176795e-17,  0.00000000e+00,
                          0.00000000e+00],
                        [-7.85046229e-17,  7.07106781e-01,  0.00000000e+00,
                          0.00000000e+00],
                        [ 0.00000000e+00,  0.00000000e+00, -7.07106781e-01,
                         -3.25176795e-17],
                        [ 0.00000000e+00,  0.00000000e+00, -7.85046229e-17,
                          7.07106781e-01]])

[genos]

So:

• let $A, B = Re\{UU^T\}, Im\{UU^T\}$
• if $A$ is diagonal, which covers the $A = 0$ case, proceed with just eigendecomp of $B$
• similarly if $B$ is diagonal, proceed with eigndecomp $A$
• otherwise, do the decomposition with both

import numpy as np
from rich import print
from rich.table import Table

def is_diag(m: np.ndarray) -> bool:
    return np.abs(m - np.diag(np.diag(m))).max() < 1e-8

def decomp(uut: np.ndarray) -> np.ndarray:
    a, b = uut.real, uut.imag
    if is_diag(a):
        _, x = np.linalg.eig(b)
    elif is_diag(b):
        _, x = np.linalg.eig(a)
    else:
        _, g = np.linalg.eig(b)
        g_inv = np.linalg.pinv(g)
        c = g_inv @ a @ g_inv.T
        _, v = np.linalg.eig(c)
        x = g_inv.T @ v.T
    return x

table = Table(title="UU^T Whack-a-Mole")
table.add_column("Unitary")
table.add_column("X^TAX Diagonal?")
table.add_column("X^TBX Diagonal?")
for name, unitary in [
    ("cnot", np.array([[1, 0, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0], [0, 1, 0, 0]])),
    ("cphase", np.diag([1, 1, 1, 1j])),
    ("latest", np.kron(np.eye(2), 1 / np.sqrt(2) * np.array([[1, -1j], [-1j, 1]]))),
]:
    a, b = unitary.real, unitary.imag
    x = decomp(unitary)
    table.add_row(name, str(is_diag(x.T @ a @ x)), str(is_diag(x.T @ b @ x)))
print(table)

               UU^T Whack-a-Mole
┏━━━━━━━━━┳━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━┓
┃ Unitary ┃ X^TAX Diagonal? ┃ X^TBX Diagonal? ┃
┡━━━━━━━━━╇━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━┩
│ cnot    │ True            │ True            │
│ cphase  │ True            │ True            │
│ latest  │ True            │ True            │
└─────────┴─────────────────┴─────────────────┘

[genos]

Nope, that just fails another test. In trying

diff --git a/src/compilers/approx.lisp b/src/compilers/approx.lisp
index 2913110..6c21563 100644
--- a/src/compilers/approx.lisp
+++ b/src/compilers/approx.lisp
@@ -179,12 +179,6 @@
                  (not (double~ 0.0d0 (magicl:tref m i j))))
         (return-from diagonal-matrix-p nil)))))

-(defun zero-matrix-p (m)
-  (dotimes (i (magicl:nrows m) t)
-    (dotimes (j (magicl:ncols m))
-      (when (not (double~ 0.0d0 (abs (magicl:tref m i j))))
-        (return-from zero-matrix-p nil)))))
-
 (defun real->complex (m)
   "Convert a real matrix M to a complex one."
   (let ((cm (magicl:zeros
@@ -204,10 +198,10 @@ are diagonal. Return (VALUES X UU^T).
          (a (magicl:map #'realpart uut))
          (b (magicl:map #'imagpart uut)))
     (cond
-      ((zero-matrix-p a)
+      ((diagonal-matrix-p a)
        (values (nth-value 1 (magicl:eig b))
                uut))
-      ((zero-matrix-p b)
+      ((diagonal-matrix-p b)
        (values (nth-value 1 (magicl:eig a))
                uut))
       (t

we fail

CL-QUIL-TESTS (Suite)
  TEST-LOGICAL-MATRIX-SANITY
    I 0
    X 0
    Y 0
    Z 0
    X 0;X 0;
    I 0;I 1;
    CNOT 0 1
    CNOT 1 0
    X 0;CNOT 0 1;
    X 0;CNOT 1 0;
    PHASE(pi/2) 0
    PHASE(-pi/2) 0
    DAGGER PHASE(pi/2) 0
    DAGGER DAGGER PHASE(pi/2) 0
    CONTROLLED X 0 1
    CONTROLLED X 1 0
    CONTROLLED Y 0 1
    CONTROLLED Y 1 0
    CONTROLLED DAGGER PHASE(pi/2) 1 0
    DAGGER CONTROLLED PHASE(pi/2) 1 0
    FORKED Y 0 1
    FORKED Y 1 0
    FORKED PHASE(pi, pi/2) 1 0
    FORKED PHASE(pi/2, pi) 1 0

with

  X^T (UU^T) X not diagonal!

X =
#<MATRIX/COMPLEX-DOUBLE-FLOAT (4x4):
  -0.942 + 0.000j     0.327 + 0.000j    -0.068 + 0.000j     0.040 + 0.000j
   0.329 + 0.000j     0.943 + 0.000j     0.003 + 0.000j     0.053 + 0.000j
  -0.005 + 0.000j    -0.052 + 0.000j     0.358 + 0.000j     0.932 + 0.000j
  -0.068 + 0.000j     0.040 + 0.000j     0.931 + 0.000j    -0.356 + 0.000j>

U =
#<MATRIX/COMPLEX-DOUBLE-FLOAT (4x4):
  -0.015 + 0.634j    -0.076 - 0.304j    -0.070 - 0.696j     0.081 + 0.053j
   0.231 - 0.046j    -0.659 - 0.097j     0.533 + 0.063j     0.459 - 0.044j
   0.304 - 0.076j     0.634 + 0.015j     0.053 - 0.081j     0.696 - 0.070j
   0.097 - 0.659j    -0.046 - 0.231j    -0.044 - 0.459j    -0.063 + 0.533j>

UU^T =
#<MATRIX/COMPLEX-DOUBLE-FLOAT (4x4):
  -0.965 + 0.133j     0.093 + 0.000j    -0.000 + 0.000j    -0.000 + 0.206j
   0.093 + 0.000j     0.965 + 0.133j     0.000 - 0.206j    -0.000 + 0.000j
  -0.000 + 0.000j     0.000 - 0.206j     0.965 - 0.133j    -0.093 + 0.000j
  -0.000 + 0.206j    -0.000 + 0.000j    -0.093 + 0.000j    -0.965 - 0.133j>

X^T(UU^T)X =
#<MATRIX/COMPLEX-DOUBLE-FLOAT (4x4):
  -0.814 + 0.159j     0.526 - 0.007j    -0.001 - 0.187j     0.028 - 0.000j
   0.526 - 0.007j     0.813 + 0.158j    -0.032 - 0.012j     0.003 - 0.187j
  -0.001 - 0.187j    -0.032 - 0.012j    -0.780 - 0.159j     0.575 + 0.007j
   0.028 - 0.000j     0.003 - 0.187j     0.575 + 0.007j     0.780 - 0.158j>

Original M such that U = E^T M E is
#<MATRIX/COMPLEX-DOUBLE-FLOAT (4x4):
  -0.466 + 0.422j    -0.019 + 0.323j    -0.144 - 0.189j     0.497 + 0.443j
   0.288 - 0.147j    -0.199 + 0.597j    -0.206 - 0.633j    -0.115 - 0.208j
   0.422 + 0.466j     0.323 + 0.019j     0.189 - 0.144j    -0.443 + 0.497j
   0.147 + 0.288j    -0.597 - 0.199j    -0.633 + 0.206j    -0.208 + 0.115j>

[stylewarning]

I think whack a mole is the only way for me, a lowly software engineer, to figure this out. (:

I'm just going to call it experimental mathematics to save face. :)

[stylewarning]

Is there a general problem for unitary matrices of the form $SU(2)\otimes SU(2)$ (which is the vast minority of $SU(4)$ )? I haven't tested yet.

[genos]

Oh, we could have $A$ or $B$ not be diagonal but skew diagonal or some other permutation matrix.

[stylewarning]

Oh, we could have $A$ or $B$ not be diagonal but skew diagonal or some other permutation matrix.

I did write a function to detect whether something looks like $\pi\cdot\mathrm{diag}(a,b,c,d)$ for a permutation $\pi$, I could dig that back up.

[stylewarning]

Pushed that here: https://github.com/quil-lang/quilc/pull/850/commits/5473335feab7c15f258a4656c9812a9a7b01cf12

[genos]

Just found out about Takagi's Factorization, which I think is what we're looking for. Math Overflow has a Python version, as does Strawberry Fields.

[genos]

import numpy as np
import scipy.linalg as la
from rich import print
from rich.table import Table

latest = np.array(
    [
        [-0.965 + 0.133j, 0.093 + 0.000j, -0.000 + 0.000j, -0.000 + 0.206j],
        [0.093 + 0.000j, 0.965 + 0.133j, 0.000 - 0.206j, -0.000 + 0.000j],
        [-0.000 + 0.000j, 0.000 - 0.206j, 0.965 - 0.133j, -0.093 + 0.000j],
        [-0.000 + 0.206j, -0.000 + 0.000j, -0.093 + 0.000j, -0.965 - 0.133j],
    ]
)

def decomp_via_takagi(m: np.ndarray) -> np.ndarray:
    """https://math.stackexchange.com/a/4448242"""
    n = m.shape[0]
    a, b = m.real, m.imag
    d, p = la.schur(np.block([[-a, b], [b, a]]))
    pos_eigenval_positions = np.diag(d) > 0
    u = p[n:, pos_eigenval_positions] + 1j * p[:n, pos_eigenval_positions]
    return la.pinv(u).T

def is_diag(m: np.ndarray) -> bool:
    return np.abs(m - np.diag(np.diag(m))).max() < 1e-8

table = Table(title="UU^T Whack-a-Mole")
table.add_column("Unitary")
table.add_column("X^T UU^T X Diagonal?")
for name, unitary in [
    ("cnot", np.array([[1, 0, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0], [0, 1, 0, 0]])),
    ("cphase", np.diag([1, 1, 1, 1j])),
    ("kron", np.kron(np.eye(2), 1 / np.sqrt(2) * np.array([[1, -1j], [-1j, 1]]))),
    ("latest", latest),
]:
    a, b = unitary.real, unitary.imag
    x = decomp_via_takagi(unitary)
    table.add_row(name, str(is_diag(x.T @ unitary @ x)))
print(table)

 python decomp.py
        UU^T Whack-a-Mole
┏━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━┓
┃ Unitary ┃ X^T UU^T X Diagonal? ┃
┡━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━┩
│ cnot    │ True                 │
│ cphase  │ True                 │
│ kron    │ True                 │
│ latest  │ True                 │
└─────────┴──────────────────────┘

[stylewarning]

@genos good sleuthing! looking into it

[genos]

Here's hoping some randomized testing will give us confidence…

import numpy as np
from rich import print
from rich.progress import track
import scipy.linalg as la
from scipy.stats import unitary_group
import typer

def decomp_via_takagi(m: np.ndarray) -> np.ndarray:
    n = m.shape[0]
    a, b = m.real, m.imag
    d, p = la.schur(np.block([[-a, b], [b, a]]))
    pos_eigenval_positions = np.diag(d) > 0
    u = p[n:, pos_eigenval_positions] + 1j * p[:n, pos_eigenval_positions]
    return la.pinv(u).T

def assert_is_diag(m: np.ndarray):
    assert (d := np.abs(m - np.diag(np.diag(m))).max()) < 1e-8, f"Max-abs off-diag: {d}"

def test_su4(iterations: int, seed: int):
    rng = unitary_group(dim=4, seed=seed)
    for _ in track(range(iterations), description="[cyan]Testing SU(4)…[/cyan]"):
        u = rng.rvs()
        uut = u @ u.T
        x = decomp_via_takagi(uut)
        assert_is_diag(x.T @ uut @ x)

def test_su2_x_su2(iterations: int, seed: int):
    rng = unitary_group(dim=2, seed=seed)
    for _ in track(range(iterations), description="[cyan]Testing SU(2)⊗SU(2)…[/cyan]"):
        u = np.kron(rng.rvs(), rng.rvs())
        uut = u @ u.T
        x = decomp_via_takagi(uut)
        assert_is_diag(x.T @ uut @ x)

def main(
    iterations: int = typer.Option(100_000, help="Number of tests to run"),
    su4_seed: int = typer.Option(96692877, help="Seed for SU(4) tests (random.org)"),
    su2_seed: int = typer.Option(29676226, help="Seed for SU(2) tests (random.org)"),
):
    test_su4(iterations, su4_seed)
    test_su2_x_su2(iterations, su2_seed)
    print("[bold][green]Passed![/green][/bold]")

if __name__ == "__main__":
    typer.run(main)

[stylewarning]

@genos just added magicl:schur in https://github.com/quil-lang/magicl/pull/182

Will try to add the above algo and see how it goes.

[stylewarning]

@genos, the line

    pos_eigenval_positions = np.diag(d) > 0
    u = p[n:, pos_eigenval_positions] + 1j * p[:n, pos_eigenval_positions]

Is this saying that if the $k$ th eigenvalue is positive (presumably there will be $n$ of them), then p[n:, ...] is an $n\times n$ matrix whose rows are the last $n$ rows of $p$, and whose columns are the positions (e.g., $k$) of positive eigenvalues?

Are we guaranteed the number of positive eigenvalues is $n$?

[stylewarning]

@genos So, as is sometimes typical with me, I started implementing things without thinking about it too deeply first. Only after quilc tests started failing did I think about it. :S

Don't we have a problem with Takagi? It just says that if we have a symmetric matrix $A$, then it will produce a unitary $V$ and a diagonal $D$ such that $VDV^{T} = A$. But if $A = UU^{T}$, then just set $V := U$ and $D := I$. Indeed, we see this immediately:

CL-QUIL> (let ((u (random-unitary 4)))
           (multiple-value-bind (v uut)
               (takagi-decomposition-of-uu^t u)
             (magicl:@ (magicl:transpose v) uut v)))
#<MAGICL:MATRIX/COMPLEX-DOUBLE-FLOAT (4x4):
   1.000 - 0.000j     0.000 + 0.000j    -0.000 + 0.000j     0.000 - 0.000j
   0.000 + 0.000j     1.000 - 0.000j    -0.000 - 0.000j    -0.000 - 0.000j
  -0.000 + 0.000j    -0.000 - 0.000j     1.000 - 0.000j    -0.000 + 0.000j
   0.000 - 0.000j    -0.000 - 0.000j     0.000 + 0.000j     1.000 - 0.000j>

We need a $UU^{T} = ODO^{T}$ with special orthogonal $O$.

[genos]

@stylewarning

1. I think we get $n$ positive eigenvalues by construction. The Math Overflow code I stole seems to set it up that way, but I haven't worked through it in detail :-/
2. Oh no! Well of course $A = UU^T \implies V \leftarrow U \;\mathrm{ and }\;D \leftarrow I$ would work 🤦. Back to internet sleuthing.

[genos]

Wait, if we're looking for a decomposition of a symmetric unitary $UU^T$, won't the eigen-decomposition suffice?

import numpy as np
from rich import print
from rich.progress import track
from scipy.stats import unitary_group
import typer

def decomp_via_eigendecomp(uut: np.ndarray) -> np.ndarray:
    return np.linalg.eig(uut)[1]

def _test(u: np.ndarray):
    uut = u @ u.T
    x = decomp_via_eigendecomp(uut)
    m = x.T @ uut @ x
    assert (d := np.abs(m - np.diag(np.diag(m))).max()) < 1e-8, f"Max-abs off-diag: {d}"

def test_su4(iterations: int, seed: int):
    rng = unitary_group(dim=4, seed=seed)
    for _ in track(range(iterations), description="[cyan]Testing SU(4)…[/cyan]"):
        _test(rng.rvs())

def test_su2_x_su2(iterations: int, seed: int):
    rng = unitary_group(dim=2, seed=seed)
    for _ in track(range(iterations), description="[cyan]Testing SU(2)⊗SU(2)…[/cyan]"):
        _test(np.kron(rng.rvs(), rng.rvs()))

def main(
    iterations: int = typer.Option(100_000, help="Number of tests to run"),
    su4_seed: int = typer.Option(96692877, help="Seed for SU(4) tests (random.org)"),
    su2_seed: int = typer.Option(29676226, help="Seed for SU(2) tests (random.org)"),
):
    test_su4(iterations, su4_seed)
    test_su2_x_su2(iterations, su2_seed)
    print("[bold][green]Passed![/green][/bold]")

if __name__ == "__main__":
    typer.run(main)

~/tmp ∃ python test_decomp.py
Testing SU(4)… ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 100% 0:00:00
Testing SU(2)⊗SU(2)… ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 100% 0:00:00
Passed!

[genos]

No, of course it's not that easy

  The assertion (CL-QUIL.FRONTEND:DOUBLE= 1.0d0 (MAGICL:DET QUIL::EVECS))
  failed with (MAGICL:DET QUIL::EVECS) =
  #C(0.9999925406226265d0 0.0038624731849657354d0).

[genos]

@stylewarning I'm having trouble understanding why the following doesn't set the determinant of the returned $X$ a.k.a QUIL::EVECS to 1:

diff --git a/src/compilers/approx.lisp b/src/compilers/approx.lisp
index 53abb86..de30a87 100644
--- a/src/compilers/approx.lisp
+++ b/src/compilers/approx.lisp
@@ -212,39 +212,13 @@
     (magicl::map-to #'complex m cm)
     cm))

-(defun takagi-decomposition-of-uu^t (u)
-  "Given a unitary U, finds an X such that
-
-    X^T (UU^T) X
-
-is a diagonal matrix. Return (VALUES X UU^T)."
+(defun decomposition-of-uu^t (u)
+  "Given a unitary U, finds an X such that X^T (UU^T) X is a diagonal matrix. Return (VALUES X UU^T)."
   (let* ((uut (magicl:@ u (magicl:transpose u)))
-         (n (magicl:nrows uut))
-         (a (magicl:.realpart uut))
-         (b (magicl:.imagpart uut))
-         (m (magicl:block-matrix (list (magicl:map #'- a) b b a) '(2 2))))
-    (multiple-value-bind (p d) (magicl:schur m)
-      (let* ((positive-eigs (loop :for i :below (magicl:nrows d)
-                                  :for e := (magicl:tref d i i)
-                                  :when (and (double~ 0.0d0 (imagpart e))
-                                             (plusp (realpart e)))
-                                    :collect i))
-             (diagonalizer (magicl:zeros (list n n) :type '(complex double-float))))
-        (assert (= 4 (length positive-eigs))
-            ()
-            "Expected 4 positive eigenvalues. Got ~D." (length positive-eigs))
-        (dotimes (row n)
-          (loop :for col :below n
-                :for from-col :in positive-eigs :do
-                  (setf (magicl:tref diagonalizer row col)
-                        (complex (magicl:tref p (+ n row) from-col)
-                                 (magicl:tref p row from-col)))))
-        (setf diagonalizer (magicl:transpose (magicl:inv diagonalizer)))
-        (when *check-math*
-          (assert (diagonal-matrix-p (magicl:@ (magicl:transpose diagonalizer)
-                                               uut
-                                               diagonalizer))))
-        (values diagonalizer uut)))))
+         (x (nth-value 1 (magicl:eig uut)))
+         (det-x (magicl:det x))
+         (normalized-x (magicl:map #'(lambda (v) (/ v det-x)) x)))
+    (values normalized-x uut)))

[stylewarning]

consider for example the diagonal matrix $M = \mathrm{diag}(\lambda_1, \ldots, \lambda_n)$, then $\det M = \lambda_1\cdots\lambda_n$, so simply dividing each entry is going to lead to $\det (M/k) = \frac{1}{k^n}\lambda_1\cdots\lambda_n$. I think you'd need $(\det x)^{-\dim x}$, if your sole goal was to normalize the determinant.

Suppose we have a different matrix $O$ which is presumed $OO^T = I$. Then $(O/k)(O/k)^{T} \neq I$.

We need both that $\det O = 1$ and $OO^T=I$.

[genos]

Ugh, the (many) professors who tried to teach me (numerical) linear algebra are very disappointed in me

[stylewarning]

@genos Maybe it's good to take a step back and re-evaluate what we might be solving by avoiding calculating eigenvalues in the current way. We've at least removed all non-determinism. I think the only benefit could be that we might get more numerical stability, but even for that I don't have a good argument.

There are other things which show that we have some misunderstandings in the code, which may be good sources for fixing. For instance, when testing, we get:

WARNING:
   Complex determinant found for a matrix expected to be (real) orthogonal: det=#C(-0.1572270927935247d0 -0.9875624746271482d0)

which is something the code is expecting, yet we don't see follow-up failures from this violation.

[stylewarning]

Some stats:

Running the test suite in full, it takes find-diagonalizer-in-e-basis

• 1 attempt to find 121,000 decompositions
• 2 attempts to find 50 decompositions
• 3+ attempts for no decompositions

So it seems the current state of affairs might be OK, at least in terms of it doing what it's supposed to.

[genos]

Sorry for going at this in a rather headstrong fashion, I was under a deadline and wanted to try to squeeze out a result beforehand.

Glad you’re taking a fresh look at things @stylewarning! When you say “the current state of affairs,” is that the PR branch or master? What ultimately brought me here was master failing to compile on ARM.

[stylewarning]

@genos ah yeah, I forgot about ARM completely. What's the state of ARM on master?

[genos]

I haven’t fuzzed the random state of the programs I was trying to generate, compile, and run, but #842 was uniformly throwing violent errors my way.

[stylewarning]

@genos If I had an ARM machine, I would be happy to do some spelunking. Unfortunately I don't. :/

Actually, the last error you posted in the thread is somewhat hopeful, it's talking about a diagonal matrix of $1/\sqrt{2}$ which seems "tame" enough to find the root of the issue.

If I had to put money on a particular thing being a problem, I'd say it might be a different in how eig behaves on each platform (?).

[genos]

@stylewarning I’d charged ahead on this so hard, I hadn’t focused on the thing that brought me here in the first place! Though I’d be embarrassed if it turned out to be that simple a fix, if all we needed was some special case handling around “is this already diagonal?” I’d be much obliged (and humbled)