`optimize_for_target_gateset` fails after `merge_k_qubit_unitaries` with ValueError: "Coefficient is not unitary"

[YaegakiY]

Describe the issue

When using merge_k_qubit_unitaries followed by optimize_for_target_gateset, the simulator raises a ValueError: Coefficient is not unitary.
The issue disappears if expand_composite is applied before merge_k_qubit_unitaries.

Explain how to reproduce the bug or problem

To reproduce:

import cirq
from cirq.transformers import *

q = cirq.LineQubit.range(6)
circuit = cirq.Circuit()


circuit.append(cirq.H(q[1]).controlled_by(q[2]))


circuit.append(cirq.measure(q, key="m"))


# circuit = expand_composite(circuit)
circuit = merge_k_qubit_unitaries(circuit, k=2)
circuit = optimize_for_target_gateset(circuit)

simulator = cirq.Simulator() 
result = simulator.run(circuit, repetitions=1)
print(result.histogram(key='m'))

Traceback (most recent call last): File "C:\Users\temp_test.py", line 16, in circuit = optimize_for_target_gateset(circuit) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "C:\Users\venv\Lib\site-packages\cirq\transformers\transformer_api.py", line 367, in func_with_logging return _transform_and_log( ^^^^^^^^^^^^^^^^^^^ File "C:\Users\venv\Lib\site-packages\cirq\transformers\transformer_api.py", line 426, in _transform_and_log transformed_circuit = _run_transformer_on_circuit( ^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "C:\Users\venv\Lib\site-packages\cirq\transformers\transformer_api.py", line 412, in _run_transformer_on_circuit return func(mutable_circuit if mutable_circuit else circuit, **kwargs) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "C:\Users\venv\Lib\site-packages\cirq\transformers\optimize_for_target_gateset.py", line 138, in optimize_for_target_gateset return _decompose_operations_to_target_gateset( ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "C:\Users\venv\Lib\site-packages\cirq\transformers\transformer_api.py", line 367, in func_with_logging return _transform_and_log( ^^^^^^^^^^^^^^^^^^^ File "C:\Users\venv\Lib\site-packages\cirq\transformers\transformer_api.py", line 426, in _transform_and_log transformed_circuit = _run_transformer_on_circuit( ^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "C:\Users\venv\Lib\site-packages\cirq\transformers\transformer_api.py", line 412, in _run_transformer_on_circuit return func(mutable_circuit if mutable_circuit else circuit, **kwargs) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "C:\Users\venv\Lib\site-packages\cirq\transformers\optimize_for_target_gateset.py", line 92, in _decompose_operations_to_target_gateset return transformer_primitives.map_operations_and_unroll( ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "C:\Users\venv\Lib\site-packages\cirq\transformers\transformer_primitives.py", line 273, in map_operations_and_unroll return _map_operations_impl( ^^^^^^^^^^^^^^^^^^^^^ File "C:\Users\venv\Lib\site-packages\cirq\transformers\transformer_primitives.py", line 192, in _map_operations_impl mapped_ops = apply_map_func(op, idx) ^^^^^^^^^^^^^^^^^^^^^^^ File "C:\Users\venv\Lib\site-packages\cirq\transformers\transformer_primitives.py", line 164, in apply_map_func mapped_ops = [*ops.flatten_to_ops(map_func(op, idx))] ^^^^^^^^^^^^^^^^^ File "C:\Users\venv\Lib\site-packages\cirq\transformers\optimize_for_target_gateset.py", line 81, in map_func return dp.decompose( ^^^^^^^^^^^^^ File "C:\Users\venv\Lib\site-packages\cirq\protocols\decompose_protocol.py", line 315, in decompose return [*_decompose_dfs(val, args)] ^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "C:\Users\venv\Lib\site-packages\cirq\protocols\decompose_protocol.py", line 207, in _decompose_dfs decomposed = decompose_once(item, default=None, flatten=False, context=args.context) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "C:\Users\venv\Lib\site-packages\cirq\protocols\decompose_protocol.py", line 375, in decompose_once decomposed = NotImplemented if method is None else method(*args, **kwargs, context=context) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "C:\Users\venv\Lib\site-packages\cirq\ops\gate_operation.py", line 157, in _decompose_with_context_ return protocols.decompose_once_with_qubits( ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "C:\Users\venv\Lib\site-packages\cirq\protocols\decompose_protocol.py", line 456, in decompose_once_with_qubits return decompose_once(val, default, tuple(qubits), flatten=flatten, context=context) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "C:\Users\venv\Lib\site-packages\cirq\protocols\decompose_protocol.py", line 378, in decompose_once decomposed = NotImplemented if method is None else method(*args, **kwargs) ^^^^^^^^^^^^^^^^^^^^^^^ File "C:\Users\venv\Lib\site-packages\cirq\ops\matrix_gates.py", line 179, in _decompose_ decomposed.append(global_phase_op.global_phase_operation(phase_delta)) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "C:\Users\venv\Lib\site-packages\cirq\ops\global_phase_op.py", line 131, in global_phase_operation return GlobalPhaseGate(coefficient, atol)() ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "C:\Users\venv\Lib\site-packages\cirq\ops\global_phase_op.py", line 35, in __init__ raise ValueError(f'Coefficient is not unitary: {coefficient!r}') ValueError: Coefficient is not unitary: np.complex128(0.5000000000000002-0.5000000000000001j)

Tell us the version of Cirq where this happens

1.6.0.dev20250702012506

[daxfohl]

I'm unable to repro this locally on today's main, on 1.6, on 1.5, or in the colab https://quantumai.google/cirq/start/start

[YaegakiY]

Thanks for looking into this!
We tested the same script on multiple environments and found that the bug occurs non-deterministically, possibly depending on the platform.
We reproduce the bug on:

• Apple M1 Pro, macOS Sequoia 15.5, cirq version: 1.6.0.dev20250702012506 1.6.0.dev20250702211226 1.6.0.dev20250704055727 1.6.0.dev20250702182429
• Apple M2 Pro, macOS Sonoma Version 14.4 (23E214), cirq version: 1.6.0.dev20250702012506 1.6.0.dev20250704055727
• Intel i7-12700KF (x86_64), Windows 11 Enterprise — 24H2 (Build 26100.4351), cirq version: 1.6.0.dev20250702012506 1.6.0.dev20250704055727

Please let us know if you'd like a full pip freeze from any of these setups.

[daxfohl]

Weird. I only have ubuntu, so someone else will have to investigate.

IDK if a pip freeze would be useful, but it sounds like you've got a good automated system going. If you could set up a git bisect job to find the offending change, then that would be very useful. From the stack trace, it seems to be something to do with matrix gate decomposition to global phase, which has had some recent work. #7118, #7283, or #7405 seem most likely.

My best guess is that the error is in #7118's implementation of phase_delta is incorrect on those platforms (though no idea why it would be), and then one of the other two PRs caused this function to get used in this particular transformation. But that's just a guess.

[YaegakiY]

Heya!

We manually bisected the issue and found that it first appears in version 1.6.0.dev20250625152617.
Earlier versions do not raise the error.

Based on the release timestamp, this seems to coincide with #7383, which modifies how global phase is handled during decomposition of controlled gates. This change is likely related to the observed behavior.

[daxfohl]

Interesting, are the circuits produced by the merge_k_qubit_unitaries transformer different across OSes, or is the only divergence in the behavior of optimize_for_target_gateset?

[YaegakiY]

We discovered that the phase_delta function in cirq/linalg/transformations.py produces incorrect results on certain platforms, and the root cause is the floating-point precision error during computation.

The test program is as follows:

import numpy as np

u1 = np.array([
    [1, 0, 0, 0],
    [0, 1, 0, 0],
    [0, 0, 1/np.sqrt(2), 1/np.sqrt(2)],
    [0, 0, 1/np.sqrt(2), -1/np.sqrt(2)],
], dtype=np.complex128)

phase = np.exp(1j * np.pi / 4)  # e^{iπ/4}
u2 = u1 * phase

def phase_delta(u1, u2):
    max_index = np.unravel_index(np.abs(u1).argmax(), u1.shape)
    return u2[max_index] / u1[max_index]

print("abs(phase_delta):", np.abs(phase_delta(u1, u2)))
print("Computed phase_delta:", phase_delta(u1, u2))

On correct platforms, the output is:

abs(phase_delta): 1.0
Computed phase_delta: (0.7071067811865476+0.7071067811865475j)

On incorrect platforms, the output is:

abs(phase_delta): 1.0000000000000002
Computed phase_delta: (0.7071067811865476+0.7071067811865476j)

Although the difference is subtle, we suspect it affects the final behavior of phase_delta.

To further investigate, we added debugging to the function:

def phase_delta(u1: np.ndarray, u2: np.ndarray) -> complex:
    """Calculates the phase delta of two unitaries.

    The delta is from u1 to u2. i.e. u1 * phase_delta(u1, u2) == u2.

    Assumes but does not verify that inputs are valid unitaries and differ only
    by phase.
    """
    # All cells will have the same phase difference. Just choose the cell with the largest
    # absolute value, to minimize rounding error.
    max_index = np.unravel_index(np.abs(u1).argmax(), u1.shape)
    print("==Delta Check==")
    print("np.abs(u1): ", np.abs(u1))
    print("np.abs(u2): ", np.abs(u2))
    print("u1.shape: ", u1.shape)
    print("max_index: ", max_index)
    print("delta: ", u2[max_index] / u1[max_index])
    print("compare:", np.abs(u1)[0,0] < np.abs(u1)[2,2])
    return u2[max_index] / u1[max_index]

The output on the failing platform is:

==Delta Check==
np.abs(u1):  [[1.00000000e+00 0.00000000e+00 3.65067599e-16 0.00000000e+00]
 [0.00000000e+00 7.07106781e-01 0.00000000e+00 7.07106781e-01]
 [4.47545209e-16 0.00000000e+00 1.00000000e+00 0.00000000e+00]
 [0.00000000e+00 7.07106781e-01 0.00000000e+00 7.07106781e-01]]
np.abs(u2):  [[1.         0.         0.         0.        ]
 [0.         1.         0.         0.        ]
 [0.         0.         0.70710678 0.70710678]
 [0.         0.         0.70710678 0.70710678]]
u1.shape:  (4, 4)
max_index:  (np.int64(2), np.int64(2))
delta:  (0.5000000000000002-0.5000000000000001j)
compare: True

The output on the correct platform is:

==Delta Check==
np.abs(u1):  [[1.00000000e+00 0.00000000e+00 2.55893763e-16 0.00000000e+00]
 [0.00000000e+00 7.07106781e-01 0.00000000e+00 7.07106781e-01]
 [2.55893763e-16 0.00000000e+00 1.00000000e+00 0.00000000e+00]
 [0.00000000e+00 7.07106781e-01 0.00000000e+00 7.07106781e-01]]
np.abs(u2):  [[1.         0.         0.         0.        ]
 [0.         1.         0.         0.        ]
 [0.         0.         0.70710678 0.70710678]
 [0.         0.         0.70710678 0.70710678]]
u1.shape:  (4, 4)
max_index:  (0, 0)
delta:  (0.7071067811865475-0.7071067811865477j)
compare: False

So on the failing platform, np.abs(u1)[0,0] < np.abs(u1)[2,2] evaluates to True, even though both are theoretically equal to 1.0. This results in (2, 2) being selected as max_index, which leads to the wrong phase_delta of (0.5000000000000002-0.5000000000000001j), instead of the correct (0.7071067811865475-0.7071067811865477j).

Conclusion:
Due to minor floating-point rounding error, the computed maximum entry changes across platforms, even though all candidates should be equal in magnitude. This incorrect choice leads to an inaccurate phase_delta, which in turn causes incorrect gate behavior downstream.

[YaegakiY]

BTW, we also tested the comparison print("compare:", np.abs(u1)[0,0] > np.abs(u1)[2,2]) and found that on the correctly behaving platforms, it returns True — even though we originally assumed the two values would be equal. There is a very subtle difference between the two values on different platforms.

[daxfohl]

Hmm, as the docstring says, the phase_delta assumes that the input unitaries differ only by phase. But it seems like here, something is calling it with two unitaries that differ by more than just phase.

The input values are coming from https://github.com/quantumlib/Cirq/blob/main/cirq-core/cirq/ops/matrix_gates.py#L174. So somehow the circuit that's created using the decomposition of the matrix has a different unitary (beyond global phase) from the matrix itself.

My guess is that the identity gate isn't preserving the qubit order as expected. I know there has been problems with circuit qubit ordering being nondeterministic in the past.