quilc
quilc copied to clipboard
quil-lang
seeking to simplify `find-diagonalizer-in-e-basis`
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
qzinto MAGICL.
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.
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.
There's also NUMERICAL METHODS FOR SIMULTANEOUS DIAGONALIZATION, which this rando Matlab file claims to implement.
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.
Since
U U^Tis symmetric, Golub and Van Loan (pg. 500 in the 4th edition) may have what we need (though I don't know ifB = 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.
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)
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)))))))
@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.
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]
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?
- Thanks @stylewarning and @kilimanjaro for working so hard on this with me!
- For cases such as the above where one or both parts are singular, is there a better approach? And how well does
MAGICLsupport checking for singularity/looking at condition numbers, etc.? - OMG GitHub markdown supports $\LaTeX$ now!?
-
Of course!
-
It's pretty scrappy. We could try to improve things here.
-
Yes! Though it's still a little rough (e.g. sometimes needing to double escape, sometimes not).
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>
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?
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
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>
Python is OK, and good catch on the CPHASE typo. My bad. I'll look into your proposed change.
I suppose we’d probably need a second special case for if the matrix is purely imaginary.
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]
In this case $UU^T = I\otimes\frac{1}{\sqrt{2}}\begin{pmatrix}1 & -i \\ -i & 1\end{pmatrix}$.
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]])
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 │
└─────────┴─────────────────┴─────────────────┘
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>
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. :)
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.
Oh, we could have $A$ or $B$ not be diagonal but skew diagonal or some other permutation matrix.
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.
Pushed that here: https://github.com/quil-lang/quilc/pull/850/commits/5473335feab7c15f258a4656c9812a9a7b01cf12