Compensate for global phase in MatrixGate decomposition (#7118)

* Change `assert_allclose` to `equal_upto_global_phase` in test, due to MatrixGate decomposition.

* Add compensating global phase op if needed to MatrixGate decomposition.

* Fix ionq test

* address PR comments

* Extract phase_delta

Author: daxfohl

cirq-core/cirq/linalg/__init__.py

@@ -82,15 +82,16 @@
 from cirq.linalg.transformations import (
     apply_matrix_to_slices as apply_matrix_to_slices,
     density_matrix_kronecker_product as density_matrix_kronecker_product,
+    can_numpy_support_shape as can_numpy_support_shape,
     match_global_phase as match_global_phase,
     partial_trace as partial_trace,
     partial_trace_of_state_vector_as_mixture as partial_trace_of_state_vector_as_mixture,
+    phase_delta as phase_delta,
     reflection_matrix_pow as reflection_matrix_pow,
     state_vector_kronecker_product as state_vector_kronecker_product,
     sub_state_vector as sub_state_vector,
     targeted_conjugate_about as targeted_conjugate_about,
     targeted_left_multiply as targeted_left_multiply,
     to_special as to_special,
     transpose_flattened_array as transpose_flattened_array,
-    can_numpy_support_shape as can_numpy_support_shape,
 )

cirq-core/cirq/linalg/transformations.py

@@ -814,3 +814,17 @@ def _can_numpy_support_dims(num_dims: int) -> bool:
 def can_numpy_support_shape(shape: Sequence[int]) -> bool:
     """Returns whether numpy supports the given shape or not numpy/numpy#5744."""
     return min(shape, default=0) >= 0 and _can_numpy_support_dims(len(shape))
+
+
+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)
+    return u2[max_index] / u1[max_index]

cirq-core/cirq/linalg/transformations_test.py

@@ -653,3 +653,11 @@ def test_transpose_flattened_array(num_dimensions):
 @pytest.mark.parametrize('shape, result', [((), True), (30 * (1,), True), ((-3, 1, -1), False)])
 def test_can_numpy_support_shape(shape: tuple[int, ...], result: bool) -> None:
     assert linalg.can_numpy_support_shape(shape) is result
+
+
+@pytest.mark.parametrize('coeff', [1, 1j, -1, -1j, 1j**0.5, 1j**0.3])
+def test_phase_delta(coeff):
+    u1 = cirq.testing.random_unitary(4)
+    u2 = u1 * coeff
+    np.testing.assert_almost_equal(linalg.phase_delta(u1, u2), coeff)
+    np.testing.assert_almost_equal(u1 * linalg.phase_delta(u1, u2), u2)

cirq-core/cirq/ops/controlled_gate.py

@@ -34,7 +34,6 @@
     controlled_operation as cop,
     diagonal_gate as dg,
     global_phase_op as gp,
-    matrix_gates,
     op_tree,
     raw_types,
 )
@@ -220,12 +219,6 @@ def _decompose_with_context_(
                     control_qid_shape=self.control_qid_shape,
                 ).on(*control_qubits)
                 return [result, controlled_phase_op]
-
-        if isinstance(self.sub_gate, matrix_gates.MatrixGate):
-            # Default decompositions of 2/3 qubit `cirq.MatrixGate` ignores global phase, which is
-            # local phase in the controlled variant and hence cannot be ignored.
-            return NotImplemented
-
         result = protocols.decompose_once_with_qubits(
             self.sub_gate,
             qubits[self.num_controls() :],

cirq-core/cirq/ops/controlled_gate_test.py

@@ -431,10 +431,23 @@ def test_controlled_gate_is_consistent(gate: cirq.Gate, should_decompose_to_targ
 @pytest.mark.parametrize(
     'gate',
     [
+        cirq.I,
+        cirq.GlobalPhaseGate(1),
+        cirq.GlobalPhaseGate(-1),
+        cirq.GlobalPhaseGate(1j),
         cirq.GlobalPhaseGate(1j**0.7),
+        cirq.Z,
         cirq.ZPowGate(exponent=1.2, global_shift=0.3),
+        cirq.CZ,
         cirq.CZPowGate(exponent=1.2, global_shift=0.3),
+        cirq.CCZ,
         cirq.CCZPowGate(exponent=1.2, global_shift=0.3),
+        cirq.X,
+        cirq.XPowGate(exponent=1.2, global_shift=0.3),
+        cirq.CX,
+        cirq.CXPowGate(exponent=1.2, global_shift=0.3),
+        cirq.CCX,
+        cirq.CCXPowGate(exponent=1.2, global_shift=0.3),
     ],
 )
 @pytest.mark.parametrize(
@@ -476,10 +489,9 @@ def _test_controlled_gate_is_consistent(
         shape = cirq.qid_shape(cgate)
         qids = cirq.LineQid.for_qid_shape(shape)
         decomposed = cirq.decompose(cgate.on(*qids))
-        if len(decomposed) < 1000:  # CCCCCZ rounding error explodes
-            first_op = cirq.IdentityGate(qid_shape=shape).on(*qids)  # To ensure same qid order
-            circuit = cirq.Circuit(first_op, *decomposed)
-            np.testing.assert_allclose(cirq.unitary(cgate), cirq.unitary(circuit), atol=1e-1)
+        first_op = cirq.IdentityGate(qid_shape=shape).on(*qids)  # To ensure same qid order
+        circuit = cirq.Circuit(first_op, *decomposed)
+        np.testing.assert_allclose(cirq.unitary(cgate), cirq.unitary(circuit), atol=1e-13)
 
 
 def test_pow_inverse():

cirq-core/cirq/ops/matrix_gates.py

@@ -14,13 +14,13 @@
 
 """Quantum gates defined by a matrix."""
 
-from typing import Any, Dict, Iterable, Optional, Tuple, TYPE_CHECKING
+from typing import Any, Dict, Iterable, List, Optional, Tuple, TYPE_CHECKING
 
 import numpy as np
 
 from cirq import _import, linalg, protocols
 from cirq._compat import proper_repr
-from cirq.ops import phased_x_z_gate, raw_types
+from cirq.ops import global_phase_op, identity, phased_x_z_gate, raw_types
 
 if TYPE_CHECKING:
     import cirq
@@ -148,18 +148,34 @@ def _phase_by_(self, phase_turns: float, qubit_index: int) -> 'MatrixGate':
         return MatrixGate(matrix=result.reshape(self._matrix.shape), qid_shape=self._qid_shape)
 
     def _decompose_(self, qubits: Tuple['cirq.Qid', ...]) -> 'cirq.OP_TREE':
+        from cirq.circuits import Circuit
+
+        decomposed: List['cirq.Operation'] = NotImplemented
         if self._qid_shape == (2,):
-            return [
+            decomposed = [
                 g.on(qubits[0])
                 for g in single_qubit_decompositions.single_qubit_matrix_to_gates(self._matrix)
             ]
         if self._qid_shape == (2,) * 2:
-            return two_qubit_to_cz.two_qubit_matrix_to_cz_operations(
+            decomposed = two_qubit_to_cz.two_qubit_matrix_to_cz_operations(
                 *qubits, self._matrix, allow_partial_czs=True
             )
         if self._qid_shape == (2,) * 3:
-            return three_qubit_decomposition.three_qubit_matrix_to_operations(*qubits, self._matrix)
-        return NotImplemented
+            decomposed = three_qubit_decomposition.three_qubit_matrix_to_operations(
+                *qubits, self._matrix
+            )
+        if decomposed is NotImplemented:
+            return NotImplemented
+        # The above algorithms ignore phase, but phase is important to maintain if the gate is
+        # controlled. Here, we add it back in with a global phase op.
+        ident = identity.IdentityGate(qid_shape=self._qid_shape).on(*qubits)  # Preserve qid order
+        u = protocols.unitary(Circuit(ident, *decomposed)).reshape(self._matrix.shape)
+        phase_delta = linalg.phase_delta(u, self._matrix)
+        # Phase delta is on the complex unit circle, so if real(phase_delta) >= 1, that means
+        # no phase delta. (>1 is rounding error).
+        if phase_delta.real < 1:
+            decomposed.append(global_phase_op.global_phase_operation(phase_delta))
+        return decomposed
 
     def _has_unitary_(self) -> bool:
         return True

cirq-core/cirq/transformers/analytical_decompositions/three_qubit_decomposition.py

@@ -24,7 +24,7 @@
 
 def three_qubit_matrix_to_operations(
     q0: ops.Qid, q1: ops.Qid, q2: ops.Qid, u: np.ndarray, atol: float = 1e-8
-) -> Sequence[ops.Operation]:
+) -> List[ops.Operation]:
     """Returns operations for a 3 qubit unitary.
 
     The algorithm is described in Shende et al.:

cirq-ionq/cirq_ionq/ionq_gateset.py

@@ -51,6 +51,7 @@ def __init__(self, *, atol: float = 1e-8):
             cirq.YYPowGate,
             cirq.ZZPowGate,
             cirq.MeasurementGate,
+            cirq.GlobalPhaseGate,
             unroll_circuit_op=False,
         )
         self.atol = atol