Add method to get a Clifford gate and phase global phase from unitary (#5568)

Fixes #2847
Review: @viathor

Author: maffoo

cirq-core/cirq/ops/clifford_gate.py

@@ -620,6 +620,29 @@ def from_unitary(u: np.ndarray) -> Optional['SingleQubitCliffordGate']:
             _to_clifford_tableau(x_to=x_to, z_to=z_to)
         )
 
+    @classmethod
+    def from_unitary_with_global_phase(
+        cls, u: np.ndarray
+    ) -> Optional[Tuple['SingleQubitCliffordGate', complex]]:
+        """Creates Clifford gate with given unitary, including global phase.
+
+        Args:
+            u: 2x2 unitary matrix of a Clifford gate.
+
+        Returns:
+            A tuple of a SingleQubitCliffordGate and a global phase, such that
+            the gate unitary (as given by `cirq.unitary`) times the global phase
+            is identical to the given unitary `u`; or `None` if `u` is not the
+            matrix of a single-qubit Clifford gate.
+        """
+        gate = cls.from_unitary(u)
+        if gate is None:
+            return None
+        # Find the entry with the largest magnitude in the input unitary, to find
+        # the global phase difference between the input unitary and the gate unitary.
+        k = max(np.ndindex(*u.shape), key=lambda t: abs(u[t]))
+        return gate, u[k] / protocols.unitary(gate)[k]
+
     def pauli_tuple(self, pauli: Pauli) -> Tuple[Pauli, bool]:
         """Returns a tuple of a Pauli operator and a boolean.
 
@@ -730,10 +753,7 @@ def _act_on_(
     # Single Clifford Gate decomposition is more efficient than the general Tableau decomposition.
     def _decompose_(self, qubits: Sequence['cirq.Qid']) -> 'cirq.OP_TREE':
         (qubit,) = qubits
-        if self == SingleQubitCliffordGate.H:
-            return (common_gates.H(qubit),)
-        rotations = self.decompose_rotation()
-        return tuple(r.on(qubit) ** (qt / 2) for r, qt in rotations)
+        return tuple(gate.on(qubit) for gate in self.decompose_gate())
 
     def _commutes_(
         self, other: Any, *, atol: float = 1e-8
@@ -773,10 +793,29 @@ def _unitary_(self) -> np.ndarray:
             mat = protocols.unitary(op).dot(mat)
         return mat
 
+    def decompose_gate(self) -> Sequence['cirq.Gate']:
+        """Decomposes this clifford into a series of H and pauli rotation gates.
+
+        Returns:
+            A sequence of H and pauli rotation gates which are equivalent to this
+            clifford gate if applied in order. This decomposition agrees with
+            cirq.unitary(self), including global phase.
+        """
+        if self == SingleQubitCliffordGate.H:
+            return [common_gates.H]
+        rotations = self.decompose_rotation()
+        return [r ** (qt / 2) for r, qt in rotations]
+
     def decompose_rotation(self) -> Sequence[Tuple[Pauli, int]]:
-        """Returns ((first_rotation_axis, first_rotation_quarter_turns), ...)
+        """Decomposes this clifford into a series of pauli rotations.
 
-        This is a sequence of zero, one, or two rotations."""
+        Each rotation is given as a tuple of (axis, quarter_turns),
+        where axis is a Pauli giving the axis to rotate about. The
+        result will be a sequence of zero, one, or two rotations.
+
+        Note that the combined unitary effect of these rotations may
+        differ from cirq.unitary(self) by a global phase.
+        """
         x_rot = self.pauli_tuple(pauli_gates.X)
         y_rot = self.pauli_tuple(pauli_gates.Y)
         z_rot = self.pauli_tuple(pauli_gates.Z)

cirq-core/cirq/ops/clifford_gate_test.py

@@ -530,26 +530,59 @@ def test_from_unitary(clifford_gate):
     result_gate = cirq.SingleQubitCliffordGate.from_unitary(u)
     assert result_gate == clifford_gate
 
+    result_gate2, global_phase = cirq.SingleQubitCliffordGate.from_unitary_with_global_phase(u)
+    assert result_gate2 == result_gate
+    assert np.allclose(cirq.unitary(result_gate2) * global_phase, u)
+
 
 def test_from_unitary_with_phase_shift():
     u = np.exp(0.42j) * cirq.unitary(cirq.SingleQubitCliffordGate.Y_sqrt)
     gate = cirq.SingleQubitCliffordGate.from_unitary(u)
 
     assert gate == cirq.SingleQubitCliffordGate.Y_sqrt
 
+    gate2, global_phase = cirq.SingleQubitCliffordGate.from_unitary_with_global_phase(u)
+    assert gate2 == gate
+    assert np.allclose(cirq.unitary(gate2) * global_phase, u)
+
 
 def test_from_unitary_not_clifford():
     # Not a single-qubit gate.
     u = cirq.unitary(cirq.CNOT)
     assert cirq.SingleQubitCliffordGate.from_unitary(u) is None
+    assert cirq.SingleQubitCliffordGate.from_unitary_with_global_phase(u) is None
 
     # Not an unitary matrix.
     u = 2 * cirq.unitary(cirq.X)
     assert cirq.SingleQubitCliffordGate.from_unitary(u) is None
+    assert cirq.SingleQubitCliffordGate.from_unitary_with_global_phase(u) is None
 
     # Not a Clifford gate.
     u = cirq.unitary(cirq.T)
     assert cirq.SingleQubitCliffordGate.from_unitary(u) is None
+    assert cirq.SingleQubitCliffordGate.from_unitary_with_global_phase(u) is None
+
+
+@pytest.mark.parametrize(
+    "clifford_gate",
+    (
+        cirq.SingleQubitCliffordGate.I,
+        cirq.SingleQubitCliffordGate.H,
+        cirq.SingleQubitCliffordGate.X,
+        cirq.SingleQubitCliffordGate.Y,
+        cirq.SingleQubitCliffordGate.Z,
+        cirq.SingleQubitCliffordGate.X_sqrt,
+        cirq.SingleQubitCliffordGate.Y_sqrt,
+        cirq.SingleQubitCliffordGate.Z_sqrt,
+        cirq.SingleQubitCliffordGate.X_nsqrt,
+        cirq.SingleQubitCliffordGate.Y_nsqrt,
+        cirq.SingleQubitCliffordGate.Z_nsqrt,
+    ),
+)
+def test_decompose_gate(clifford_gate):
+    gates = clifford_gate.decompose_gate()
+    u = functools.reduce(np.dot, [np.eye(2), *(cirq.unitary(gate) for gate in reversed(gates))])
+    assert np.allclose(u, cirq.unitary(clifford_gate))  # No global phase difference.
 
 
 @pytest.mark.parametrize('trans_x,trans_z', _all_rotation_pairs())

cirq-core/cirq/sim/clifford/stabilizer_simulation_state.py

@@ -161,21 +161,14 @@ def _strat_act_from_single_qubit_decompose(
             if not has_unitary(val):
                 return NotImplemented
             u = unitary(val)
-            clifford_gate = SingleQubitCliffordGate.from_unitary(u)
-            if clifford_gate is not None:
-                # Gather the effective unitary applied so as to correct for the
-                # global phase later.
-                final_unitary = np.eye(2)
-                for axis, quarter_turns in clifford_gate.decompose_rotation():
-                    gate = axis ** (quarter_turns / 2)
+            gate_and_phase = SingleQubitCliffordGate.from_unitary_with_global_phase(u)
+            if gate_and_phase is not None:
+                clifford_gate, global_phase = gate_and_phase
+                # Apply gates.
+                for gate in clifford_gate.decompose_gate():
                     self._strat_apply_gate(gate, qubits)
-                    final_unitary = np.matmul(unitary(gate), final_unitary)
-
-                # Find the entry with the largest magnitude in the input unitary.
-                k = max(np.ndindex(*u.shape), key=lambda t: abs(u[t]))
-                # Correct the global phase that wasn't conserved in the above
-                # decomposition.
-                self._state.apply_global_phase(u[k] / final_unitary[k])
+                # Apply global phase.
+                self._state.apply_global_phase(global_phase)
                 return True
 
         return NotImplemented