instability in kak decomposition

[babbush]

In some cases the kak decomposition appears to have an instability in the sign of the ZZ coefficient. Consider the following example:

a, b = cirq.LineQubit(0), cirq.LineQubit(1)

theta = 0.3
circuit_1 = cirq.Circuit()
circuit_1.append(cirq.ISWAP.on(a, b))
circuit_1.append(
    cirq.ZZPowGate(exponent=2 * theta / numpy.pi, global_shift=-0.5).on(a, b))
print(cirq.kak_decomposition(cirq.unitary(circuit_1)).interaction_coefficients)

theta += 0.001
circuit_2 = cirq.Circuit()
circuit_2.append(cirq.ISWAP.on(a, b))
circuit_2.append(
    cirq.ZZPowGate(exponent=2 * theta / numpy.pi, global_shift=-0.5).on(a, b))
print(cirq.kak_decomposition(cirq.unitary(circuit_2)).interaction_coefficients)

print(numpy.amax(numpy.absolute(cirq.unitary(circuit_1) - cirq.unitary(circuit_2))))

returns

(0.7853981633974483, 0.7853981633974483, -0.2999999999999998)
(0.7853981633974482, 0.7853981633974482, 0.30099999999999993)
0.0009999999583333466

And thus we see that despite the fact the unitaries are very similar, the kak decomposition gives a seemingly random sign on the ZZ interaction coefficient. Is this a bug?

[Strilanc]

When I run this repro script I get the following results:

(0.7853981633974483, 0.7853981633974482, 0.2999999999999998)
(0.7853981633974483, 0.7853981633974482, 0.3009999999999997)
0.0009999999583333466

With no sign instability.

This happens both in my development environment and in a freshly created virtual env. I also tried it in both python 2 and python 3.

[Strilanc]

What version of cirq are you using? What OS?

[babbush]

Mac, Mojave, Python 3.7, Cirq 0.5. The discrepancy between our builds is disturbing...

[vtomole]

I just ran this on cirq==0.5.0 and python==3.7.3 on ubuntu==18.04 with no issues; so it's probably macOS problem.

[vtomole]

@babbush This also passes in macOS 10.13, Xcode 9.4.1 : https://travis-ci.com/quantumlib/Cirq/jobs/201188271. This might be a very strange corner case, but it's worth getting to the bottom of it!

[Strilanc]

@babbush I recommend putting your laptop next to your workstation, getting this code on both of them in pycharm, and debugging through the test. Step over methods until you find the point of divergence. Then step into the method where it diverges, and step over submethods to find the more exact point of divergence. Just keep narrowing it down.

[babbush]

Just to make sure @vtomole when you say "it passes" do you mean that you get the same numbers as Craig? My code should still run, it just gives a negative sign on the .2999 number. Anyways, today I ran this on my linux workstation and see the same thing. I highly suspect the problem is with scipy or numpy's diagonalization routines so the version or cirq or the os might not be the problem. For both mac and linux I am using anaconda python version 3.7. This is one of the most common python installations for scientific programming and I am certain that many of our users use it. Can one of you try running this with anaconda 3.7 and see if you get that negative sign? Realistically, I don't foresee myself having the time to track this bug down but its quite a problem for us...

[vtomole]

@babbush

When you say "it passes" do you mean that you get the same numbers as Craig?

Yes. I haven't been able to reproduce the negative sign on .2999 on macOS, Windows, or Linux. Let me try Anaconda python==3.7 on Linux.

[vtomole]

No luck on the latest version of Anaconda:

(base) vtomole@vtomolet:~$ python test3.py 
(0.7853981633974483, 0.7853981633974482, 0.2999999999999998)
(0.7853981633974483, 0.7853981633974482, 0.3009999999999997)
0.0009999999583333466
(base) vtomole@vtomole:~$ python
Python 3.7.3 (default, Mar 27 2019, 22:11:17) 
[GCC 7.3.0] :: Anaconda, Inc. on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> exit()
(base) vtomole@vtomole:~$ pip show numpy
Name: numpy
Version: 1.16.2
Summary: NumPy is the fundamental package for array computing with Python.
Home-page: https://www.numpy.org
Author: Travis E. Oliphant et al.
Author-email: None
License: BSD
Location: /home/vagrant/anaconda3/lib/python3.7/site-packages
Requires: 
Required-by: tables, seaborn, scikit-learn, PyWavelets, pytest-doctestplus, pytest-arraydiff, patsy, pandas, numexpr, numba, mkl-random, mkl-fft, matplotlib, h5py, cirq, Bottleneck, bokeh, bkcharts, astropy
(base) vtomole@vtomolet:~$ pip show scipy
Name: scipy
Version: 1.2.1
Summary: SciPy: Scientific Library for Python
Home-page: https://www.scipy.org
Author: None
Author-email: None
License: BSD
Location: /home/vagrant/anaconda3/lib/python3.7/site-packages
Requires: 
Required-by: seaborn, scikit-learn, cirq

[babbush]

So strange. Can you try running this in an anaconda jupyter notebook and see if that changes anything?

On Mon, May 20, 2019, 09:44 Victory Omole notifications@github.com wrote:

No luck on the latest version of Anaconda:

(base) vtomole@vtomolet:~$ python test3.py (0.7853981633974483, 0.7853981633974482, 0.2999999999999998) (0.7853981633974483, 0.7853981633974482, 0.3009999999999997) 0.0009999999583333466 (base) vtomole@vtomole:~$ python Python 3.7.3 (default, Mar 27 2019, 22:11:17) [GCC 7.3.0] :: Anaconda, Inc. on linux Type "help", "copyright", "credits" or "license" for more information.

exit() (base) vtomole@vtomole:~$ pip show numpy Name: numpy Version: 1.16.2 Summary: NumPy is the fundamental package for array computing with Python. Home-page: https://www.numpy.org Author: Travis E. Oliphant et al. Author-email: None License: BSD Location: /home/vagrant/anaconda3/lib/python3.7/site-packages Requires: Required-by: tables, seaborn, scikit-learn, PyWavelets, pytest-doctestplus, pytest-arraydiff, patsy, pandas, numexpr, numba, mkl-random, mkl-fft, matplotlib, h5py, cirq, Bottleneck, bokeh, bkcharts, astropy (base) vtomole@vtomolet:~$ pip show scipy Name: scipy Version: 1.2.1 Summary: SciPy: Scientific Library for Python Home-page: https://www.scipy.org Author: None Author-email: None License: BSD Location: /home/vagrant/anaconda3/lib/python3.7/site-packages Requires: Required-by: seaborn, scikit-learn, cirq

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/quantumlib/Cirq/issues/1647?email_source=notifications&email_token=ADMVG75G4BT5L5ZNFGSXHALPWLIN3A5CNFSM4HNQXFIKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGODVZNICA#issuecomment-494064648, or mute the thread https://github.com/notifications/unsubscribe-auth/ADMVG76STZDZTBNVVS7OQQLPWLIN3ANCNFSM4HNQXFIA .

[vtomole]

Same thing: https://ibb.co/ZT2Spy5

[babbush]

This is so weird. I'm not the only one though. @ncrubin is also experiencing the same bug on both mac and linux. He has some theory about how it might be caused by some confluence of anaconda installed a certain way on intel processors. I'll let him elaborate.

[ncrubin]

I am observing the same behavior as @babbush . Running MacOS Mojave, Anaconda python3.7, Cirq 0.6.0.dev, numpy 1.16.3. Also observe the same negative sign issue on debian with the same anaconda python, cirq, and numpy specs.

Name: scipy
Version: 1.2.1
Summary: SciPy: Scientific Library for Python
Home-page: https://www.scipy.org
Author: None
Author-email: None
License: BSD
Requires: numpy
Required-by: seaborn, openfermioncirq, openfermion, cirq

$ python -m pip show numpy
Name: numpy
Version: 1.16.3
Summary: NumPy is the fundamental package for array computing with Python.
Home-page: https://www.numpy.org
Author: Travis E. Oliphant et al.
Author-email: None
License: BSD
Requires: 
Required-by: seaborn, scipy, qcelemental, pyswarm, pandas, mkl-random, mkl-fft, matplotlib, h5py, openfermioncirq, openfermion, cirq

$ python -m pip show cirq
Name: cirq
Version: 0.6.0.dev0
Summary: A framework for creating, editing, and invoking Noisy Intermediate Scale Quantum (NISQ) circuits.
Home-page: http://github.com/quantumlib/cirq
Author: The Cirq Developers
Author-email: cirq@googlegroups.com
License: Apache 2
Requires: google-api-python-client, matplotlib, networkx, numpy, protobuf, requests, sortedcontainers, scipy, sympy, typing-extensions

[vtomole]

Also observe the same negative sign issue on debian with the same anaconda python, cirq, and numpy specs.

How did you run the Debian? Was it in a Virtual Machine on your MacBook?

[ncrubin]

@vtomole : I used a separate machine running Debian.

[vtomole]

I've just reproduced this on Windows.

(base) C:\Cirq\check>python
Python 3.7.0 (default, Jun 28 2018, 08:04:48) [MSC v.1912 64 bit (AMD64)] :: Anaconda, Inc. on win32
Type "help", "copyright", "credits" or "license" for more information.
>>> import cirq
>>> import numpy
>>> a, b = cirq.LineQubit(0), cirq.LineQubit(1)
>>> theta = 0.3
>>> circuit_1 = cirq.Circuit()
>>> circuit_1.append(cirq.ISWAP.on(a, b))
>>> circuit_1.append(
...     cirq.ZZPowGate(exponent=2 * theta / numpy.pi, global_shift=-0.5).on(a, b))
>>> print(cirq.kak_decomposition(cirq.unitary(circuit_1)).interaction_coefficients)
(0.7853981633974483, 0.7853981633974483, -0.2999999999999998)

[Strilanc]

I tried on my macbook but it also doesn't reproduce the issue.

@ncrubin @vtomole are you able to narrow down the issue more specifically by stepping through with a debugger and seeing where a good machine and bad machine deviate?

[ncrubin]

Currently @vtomole is the only one with a "good machine" and a "bad machine"

[vtomole]

@Strilanc I'm currently stepping through the bad and good machine at the moment.

[babbush]

awesome @vtomole ! I'll definitely owe you some good Google schwag if you can pin this one down and implement a fix. Its currently blocking me on something important!

[vtomole]

@Strilanc @ncrubin @Strilanc Got it.

https://github.com/quantumlib/Cirq/blob/5f49f39bb05eee222fa382fe5c7f7a050de9f7b3/cirq/linalg/decompositions.py#L455 goes to https://github.com/quantumlib/Cirq/blob/5f49f39bb05eee222fa382fe5c7f7a050de9f7b3/cirq/linalg/decompositions.py#L440

This loop never runs on the implementations that give you the negative sign. This is because on the implementations that give you the negative sign, v[2] = 0.7853981633974483. For the one that doesn't it's v[2] = 0.7853981633974484.

0.7853981633974484 > np.pi / 4 is True while 0.7853981633974483 > np.pi /4 is False.

This could come from anything (hardware e.t.c). Should we just round this up a decimal place or two; or rethink the implementation of this procedure?

[vtomole]

I should also add that this is coming from

https://github.com/quantumlib/Cirq/blob/5f49f39bb05eee222fa382fe5c7f7a050de9f7b3/cirq/linalg/decompositions.py#L523 Where z = 0.7853981633974483 on negative ones and z = 0.7853981633974484 on non-negatives; so where we should round is not clear, but i prefer rounding in def canonical_shift.

[vtomole]

If you round pi / 4 down, it's 0.7853981633974483. No wonder this was hard to track.

[Strilanc]

Oh, it's just because we're on the border for the angular wrap around? Blergh. So it's not actually different, it just looks different.

Given how confusing this was, I think we should adjust the canonical range to be (-pi/4+epsilon, pi/4+epsilon]. The KakDecomposition class also needs to support approximate equality.

[Strilanc]

Okay, I've managed to find a variant of this issue that I can reproduce:

    perturb = -0.00001
    op1 = yy**perturb
    op2 = np.exp(-1j * perturb * np.pi * yy)
    assert cirq.allclose_up_to_global_phase(cirq.unitary(op1),
                                            cirq.unitary(op2),
                                            atol=1e-12)

    c1 = cirq.Circuit.from_ops(cirq.ISWAP(a, b), cirq.CZ(a, b)**0.5, op1)
    c2 = cirq.Circuit.from_ops(cirq.ISWAP(a, b), cirq.CZ(a, b)**0.5, op2)
    kak1 = cirq.kak_decomposition(cirq.unitary(c1))
    kak2 = cirq.kak_decomposition(cirq.unitary(c2))
    print(kak1.interaction_coefficients)
    print(kak2.interaction_coefficients)

Prints

(0.7853981633974483, 0.7853824554341802, -0.39269908169872414)
(0.7853981633974483, 0.7853824554341803, 0.39269908169872425)

Even though op1 is approximately equal to op2, the coefficients get negated.

[Strilanc]

Apparently these two KAK decompositions represent the same unitary (up to global phase):

KAK:
    xyz*(4/π): 1, 1, -0.5
    before: (1*π around 0.924X-0.383Y) ⊗ (0.25*π around X)
    after: (1*π around X) ⊗ (1*π around X)
KAK:
    xyz*(4/π): 1, 1, 0.5
    before: (-0.75*π around X) ⊗ (0.25*π around X)
    after: (0*π around X) ⊗ (1*π around X)

I'm not sure what to make of it, since they disagree so wildly on the Z. It may have to do specifically with the fact that the X and Y coefficients are equal to 1. I seem to remember there being something about one of the border lines of the space being redundant; maybe that's what this is.

[Strilanc]

Yes, the issue is that the canonicalization is not handling a special case where X=Y=pi/4. When this occurs, commuting an X operation across the XYZ part causes the Z part to negate without affecting the other two.

This:

is equal to this:

Up to global phase.

[babbush]

Ok, I think I see what's going on. I would still consider this behavior by itself to be buggy. It's not very helpful to have a canonicalized form which spits out different representations on different machines for the same operation...

[Strilanc]

Yes, it's a bug that it doesn't canonicalize this case.

[babbush]

The bug does not appear to end there. This is also something very strange going where the single qubit kak coefficients don't correspond to what they are supposed to be. I'm sorry for the rather long example, but below you will find pretty much the smallest example I could make that demonstrates the real problem I'm facing. I would be thrilled with even a very hacky solution to the problem I'm facing here. The code below runs but when you uncomment what should be an inconsequential line, it fails.

import cirq
import numpy

# Use kak decomposition to find single qubit rotations that match.
def find_local_equivalents(unitary1: numpy.ndarray, unitary2: numpy.ndarray):
    kak_u1 = cirq.kak_decomposition(unitary1)
    kak_u2 = cirq.kak_decomposition(unitary2)

    u_0 = kak_u1.single_qubit_operations_after[0] @ kak_u2.single_qubit_operations_after[0].conj().T
    u_1 = kak_u1.single_qubit_operations_after[1] @ kak_u2.single_qubit_operations_after[1].conj().T

    v_0 = kak_u2.single_qubit_operations_before[0].conj().T @ kak_u1.single_qubit_operations_before[0]
    v_1 = kak_u2.single_qubit_operations_before[1].conj().T @ kak_u1.single_qubit_operations_before[1]
    return v_0, v_1, u_0, u_1

# Because I want a random gate.
class RandomGate(cirq.TwoQubitGate, cirq.InterchangeableQubitsGate):

    def __init__(self, theta):
        self.theta = theta

    def _decompose_(self, qubits):
        a, b = qubits

        # Set interaction part.
        circuit = cirq.Circuit()
        angle_offset = numpy.pi / 24
        circuit.append(cirq.XXPowGate(global_shift=-0.5, exponent=0.5).on(a, b))
        circuit.append(cirq.YYPowGate(global_shift=-0.5, exponent=0.5).on(a, b))
        circuit.append(cirq.CZ(a, b) ** (-1 / 6))
        circuit.append(cirq.ZZPowGate(
            exponent=2 * (self.theta - angle_offset) / numpy.pi, global_shift=-0.5).on(a, b))

        # Get the intended circuit.
        intended_circuit = cirq.Circuit()
        intended_circuit.append(cirq.ISWAP.on(a, b))
        intended_circuit.append(cirq.ZZPowGate(
            exponent=2 * self.theta / numpy.pi, global_shift=-0.5).on(a, b))

        # Get corrections.
        b_0, b_1, a_0, a_1 = find_local_equivalents(
            cirq.unitary(intended_circuit), cirq.unitary(circuit))

        # Apply initial corrections.
        for term in cirq.single_qubit_matrix_to_gates(b_0):
            yield term.on(a)
        for term in cirq.single_qubit_matrix_to_gates(b_1):
            yield term.on(b)

        # Apply interaction part.
        yield circuit

        # Apply final corrections.
        for term in cirq.single_qubit_matrix_to_gates(a_0):
            yield term.on(a)
        for term in cirq.single_qubit_matrix_to_gates(a_1):
            yield term.on(b)

# Test.
theta = 0.1860
#theta -= 0.0001 # Uncomment this line to break things
a, b = cirq.LineQubit(0), cirq.LineQubit(1)

intended_circuit = cirq.Circuit()
intended_circuit.append(cirq.ISWAP.on(a, b))
intended_circuit.append(
    cirq.ZZPowGate(exponent=2 * theta / numpy.pi, global_shift=-0.5).on(a, b))
intended_unitary = cirq.unitary(intended_circuit)
kak_intended = cirq.kak_decomposition(intended_unitary)
print(kak_intended.interaction_coefficients)

actual_circuit = RandomGate(theta).on(a, b)
actual_unitary = cirq.unitary(actual_circuit)
kak_actual = cirq.kak_decomposition(actual_unitary)
print(kak_actual.interaction_coefficients)

cirq.testing.assert_circuits_with_terminal_measurements_are_equivalent(
    cirq.Circuit.from_ops(intended_circuit),
    cirq.Circuit.from_ops(actual_circuit), atol=1e-6)

[Strilanc]

That snippet doesn't reproduce for me.

Is there a reason you expect it to be a different issue from this canonicalization issue?

[babbush]

As I understand the canonicalization issue, we're always getting a valid KAK decomposition, but sometimes there are multiple valid KAK decompositions, and in that case the code is failing to consistently choose the same one.

However, my code above suggests something beyond this is happening. I believe my code above should always work provided the KAK decomposition is valid.

[babbush]

When you say it doesn't reproduce do you mean that uncommenting the line doesn't break things? Just set theta to random values until you see it break. Most likely you can unbreak it from that point by adding some epsilon to the theta that broke it.

[Strilanc]

You should phrase your repros in the form of a failing assertion, instead of a line that needs to be uncommented. I'm not sure what is supposed to be going wrong...?

[babbush]

The code is supposed to pass whether or not that one line is commented. But on my computers, the code only passes if that line remains commented. Sorry I made my example poorly.

[Strilanc]

Oh I'm sorry, you did phrase it as a test. It just doesn't fail on my machine.

[babbush]

well presumably that's because you aren't running it on a build that reproduces the same kak instability that @ncrubin and myself experience. I would imagine that @ncrubin would experience the same problem for the snippet and that @vtomole can reproduce it on his windows build.

[c-poole]

On my machine I do not have the sign instability but the circuit comparison test does break when the line is uncommented. If I change the offending line by adding another zero to theta -= 0.00001 it passes.

[c-poole]

If we change the line that instantiates actual_circuit to actual_circuit = cirq.Circuit.from_ops(cirq.decompose(RandomGate(theta).on(a,b))) we can get a circuit diagram from the error message. This results in the below diagrams for the failing case on my machine. What does the diagram look like on a machine where un-commenting the line doesn't cause the test to fail?

[babbush]

@vtomole can you let us know the answer to @c-poole 's question?

[vtomole]

I just ran this on my windows build

import cirq
import numpy

# Use kak decomposition to find single qubit rotations that match.
def find_local_equivalents(unitary1: numpy.ndarray, unitary2: numpy.ndarray):
    kak_u1 = cirq.kak_decomposition(unitary1)
    kak_u2 = cirq.kak_decomposition(unitary2)

    u_0 = kak_u1.single_qubit_operations_after[0] @ kak_u2.single_qubit_operations_after[0].conj().T
    u_1 = kak_u1.single_qubit_operations_after[1] @ kak_u2.single_qubit_operations_after[1].conj().T

    v_0 = kak_u2.single_qubit_operations_before[0].conj().T @ kak_u1.single_qubit_operations_before[0]
    v_1 = kak_u2.single_qubit_operations_before[1].conj().T @ kak_u1.single_qubit_operations_before[1]
    return v_0, v_1, u_0, u_1

# Because I want a random gate.
class RandomGate(cirq.TwoQubitGate, cirq.InterchangeableQubitsGate):

    def __init__(self, theta):
        self.theta = theta

    def _decompose_(self, qubits):
        a, b = qubits

        # Set interaction part.
        circuit = cirq.Circuit()
        angle_offset = numpy.pi / 24
        circuit.append(cirq.XXPowGate(global_shift=-0.5, exponent=0.5).on(a, b))
        circuit.append(cirq.YYPowGate(global_shift=-0.5, exponent=0.5).on(a, b))
        circuit.append(cirq.CZ(a, b) ** (-1 / 6))
        circuit.append(cirq.ZZPowGate(
            exponent=2 * (self.theta - angle_offset) / numpy.pi, global_shift=-0.5).on(a, b))

        # Get the intended circuit.
        intended_circuit = cirq.Circuit()
        intended_circuit.append(cirq.ISWAP.on(a, b))
        intended_circuit.append(cirq.ZZPowGate(
            exponent=2 * self.theta / numpy.pi, global_shift=-0.5).on(a, b))

        # Get corrections.
        b_0, b_1, a_0, a_1 = find_local_equivalents(
            cirq.unitary(intended_circuit), cirq.unitary(circuit))

        # Apply initial corrections.
        for term in cirq.single_qubit_matrix_to_gates(b_0):
            yield term.on(a)
        for term in cirq.single_qubit_matrix_to_gates(b_1):
            yield term.on(b)

        # Apply interaction part.
        yield circuit

        # Apply final corrections.
        for term in cirq.single_qubit_matrix_to_gates(a_0):
            yield term.on(a)
        for term in cirq.single_qubit_matrix_to_gates(a_1):
            yield term.on(b)

# Test.
theta = 0.1860
theta -= 0.0001 # Uncomment this line to break things
a, b = cirq.LineQubit(0), cirq.LineQubit(1)

intended_circuit = cirq.Circuit()
intended_circuit.append(cirq.ISWAP.on(a, b))
intended_circuit.append(
    cirq.ZZPowGate(exponent=2 * theta / numpy.pi, global_shift=-0.5).on(a, b))
intended_unitary = cirq.unitary(intended_circuit)
kak_intended = cirq.kak_decomposition(intended_unitary)
print(kak_intended.interaction_coefficients)

actual_circuit = RandomGate(theta).on(a, b)
actual_unitary = cirq.unitary(actual_circuit)
kak_actual = cirq.kak_decomposition(actual_unitary)
print(kak_actual.interaction_coefficients)

cirq.testing.assert_circuits_with_terminal_measurements_are_equivalent(
    cirq.Circuit.from_ops(intended_circuit),
    cirq.Circuit.from_ops(actual_circuit), atol=1e-6)

print(actual_circuit)
print(intended_circuit)

Result

(0.7853981633974483, 0.7853981633974482, -0.18589999999999973)
(0.7853981633974483, 0.7853981633974483, -0.18590000000000018)
<__main__.RandomGate object at 0x00000224E35B8E10>(0, 1)
0: ───iSwap───ZZ─────────
      │       │
1: ───iSwap───ZZ^0.118───

When i commented theta -= 0.0001 # Uncomment this line to break things

I got

(0.7853981633974483, 0.7853981633974483, -0.18599999999999994)
(0.7853981633974483, 0.7853981633974483, -0.18599999999999994)
<__main__.RandomGate object at 0x0000028F95B03BE0>(0, 1)
0: ───iSwap───ZZ─────────
      │       │
1: ───iSwap───ZZ^0.118───

[babbush]

Oh right, of course that's what it does. In any event, do people see why I am claiming this most recent example shows a related, but I believe, distinct, failure from the canonicalization issue? The canonicalization issue is basically occurring because, due to finite precision, there are two fairly different looking kak decompositions that are "valid" to within numerical precision (at least that is @Strilanc 's hypothesis, which seems reasonable). However, this other failure I'm pointing out shows that something about the kak decomposition is actually invalid in some cases near this instability in the sense that a sequence of gates corresponding to the kak decomposition does not reproduce the original unitary to within numerical precision.

[babbush]

And yet, this most recent problem I'm pointing out only occurs when finite precision pushes you to one side of the numerical instability, and not the other. So either something in my example code, or something in the kak decomposition code, is somehow expecting a certain convention about what to do near this instability, but I cannot figure out what it is!

[c-poole]

I agree with @babbush's conclusion that there are two separate problems that have been discovered. I've tried varying the amount theta is shifted to see if there is a clear pattern with the flipped sign on the ZZ component of the kak decomposition and when the code fails. All four conditions (flipped ZZ, code passes, flipped ZZ code fails, etc.) have occurred.

[c-poole]

I think I have a fix. I adjusted the kak_canonicalize_vector method so that this minus sign issue doesn't occur and when I did that the larger test case @babbush came up with always seems to pass. I'll put in a PR for the fix as soon as I can.

[c-poole]

To elaborate on what appeared to be happening and what my current fix does: @Strilanc correctly identified the issue as being a special case when X=pi/4 (Having Y=pi/4 doesn't appear to be necessary). When this is true, you can shift X by -pi/2 and then negate both X and Z and end up with a net operation of just negating Z. The reason @babbush ran into an issue with the equivalent circuits test is that numerical precision was sometimes the difference between correctly and incorrectly applying the local operations associated with shifting X by -pi/2 and then negating X and Z. We can resolve this issue by choosing a convention that when X=pi/4, Z>=0, which is what my PR now implements after feedback from Craig.

[dabacon]

Does #1680 fix this or is there still more work to be done here?

[Strilanc]

I asked @babbush and he said it's no longer blocking him, therefore it's good to close.

and...... clo-