Fix conjugated_by of PauliString. (#7065)

* Fix conjugated_by of PauliString.

* fix comments

* Improve assert_conjugation and associated tests.

Author: babacry

cirq-core/cirq/ops/clifford_gate.py

@@ -377,7 +377,10 @@ def __init__(self, *, _clifford_tableau: qis.CliffordTableau) -> None:
         #  ZI  [ 0  0 | 1  0  | 1 ]
         #  IZ  [ 1  0 | 1  1  | 0 ]
         # Take the third row as example: this means the ZI gate after the this gate,
-        # more precisely the conjugate transformation of ZI by this gate, becomes -ZI.
+        # more precisely the conjugate transformation of ZI by this gate, becomes -ZI:
+        #   ---(CliffordGate^-1)---ZI---CliffordGate---
+        # = unitary(CliffordGate)@unitary(ZI)@unitary(CliffordGate).conj().T
+        # = -ZI.
         # (Note the real clifford tableau has to satify the Symplectic property.
         # here is just for illustration)
         object.__setattr__(self, '_clifford_tableau', _clifford_tableau.copy())
@@ -438,7 +441,7 @@ def __repr__(self) -> str:
     def _commutes_(
         self, other: Any, *, atol: float = 1e-8
     ) -> Union[bool, NotImplementedType, None]:
-        # Note even if we assume two gates define the tabluea based on the same qubit order,
+        # Note even if we assume two gates define the tableau based on the same qubit order,
         # the following approach cannot judge it:
         # self.clifford_tableau.then(other.clifford_tableau) == other.clifford_tableau.then(
         #     self.clifford_tableau

cirq-core/cirq/ops/pauli_string.py

@@ -55,6 +55,7 @@
     pauli_gates,
     pauli_interaction_gate,
     raw_types,
+    dense_pauli_string,
 )
 
 if TYPE_CHECKING:
@@ -907,14 +908,13 @@ def dense(self, qubits: Sequence[TKey]) -> 'cirq.DensePauliString':
         Raises:
             ValueError: If the number of qubits is too small.
         """
-        from cirq.ops.dense_pauli_string import DensePauliString
 
         if not self.keys() <= set(qubits):
             raise ValueError('not self.keys() <= set(qubits)')
         # pylint: disable=too-many-function-args
         pauli_mask = [self.get(q, identity.I) for q in qubits]
         # pylint: enable=too-many-function-args
-        return DensePauliString(pauli_mask, coefficient=self.coefficient)
+        return dense_pauli_string.DensePauliString(pauli_mask, coefficient=self.coefficient)
 
     def conjugated_by(self, clifford: 'cirq.OP_TREE') -> 'PauliString':
         r"""Returns the Pauli string conjugated by a clifford operation.
@@ -976,17 +976,68 @@ def conjugated_by(self, clifford: 'cirq.OP_TREE') -> 'PauliString':
         Returns:
             The Pauli string conjugated by the given Clifford operation.
         """
-        pauli_map = dict(self._qubit_pauli_map)
-        should_negate = False
-        for op in list(op_tree.flatten_to_ops(clifford))[::-1]:
-            if pauli_map.keys().isdisjoint(set(op.qubits)):
-                continue
-            for clifford_op in _decompose_into_cliffords(op)[::-1]:
-                if pauli_map.keys().isdisjoint(set(clifford_op.qubits)):
-                    continue
-                should_negate ^= _pass_operation_over(pauli_map, clifford_op, False)
-        coef = -self._coefficient if should_negate else self.coefficient
-        return PauliString(qubit_pauli_map=pauli_map, coefficient=coef)
+
+        # Initialize the ps the same as self.
+        ps = PauliString(qubit_pauli_map=self._qubit_pauli_map, coefficient=self.coefficient)
+        all_ops = list(op_tree.flatten_to_ops(clifford))
+        all_qubits = set.union(set(self.qubits), [q for op in all_ops for q in op.qubits])
+
+        # Iteratively calculate the conjugation in reverse order of ops.
+        for op in all_ops[::-1]:
+            # To calcuate the conjugation of P (`ps`) with respect to C (`op`)
+            # Decompose P = Pc⊗R, where Pc acts on the same qubits as C, R acts on the remaining.
+            # Then the conjugation = (C^{-1}⊗I·Pc⊗R·C⊗I) = (C^{-1}·Pc·C)⊗R.
+
+            # Isolate R
+            remain: 'cirq.PauliString' = PauliString()
+            for q in all_qubits:
+                pauli = ps.get(q)
+                if pauli is not None and not q in op.qubits:
+                    remain *= pauli(q)
+
+            # Initialize the conjugation of Pc.
+            conjugated: 'cirq.DensePauliString' = (
+                dense_pauli_string.DensePauliString(pauli_mask=[identity.I for _ in op.qubits])
+                * self.coefficient
+            )
+
+            # Calculate the conjugation via CliffordGate's clifford_tableau.
+            # Note the clifford_tableau in CliffordGate represents C·P·C^-1 instead of C^-1·P·C.
+            # So we take the inverse of the tableau to match the definition of the conjugation here.
+            gate_in_clifford: 'cirq.CliffordGate'
+            if isinstance(op.gate, cirq.CliffordGate):
+                gate_in_clifford = op.gate
+            else:
+                # Convert the clifford gate to CliffordGate type.
+                gate_in_clifford = clifford_gate.CliffordGate.from_op_list([op], op.qubits)
+            tableau = gate_in_clifford.clifford_tableau.inverse()
+
+            # Calculate the conjugation by `op` via mutiplying the conjugation of each Pauli:
+            #   C^{-1}·(P_1⊗...⊗P_n)·C
+            # = C^{-1}·(P_1⊗I) ...·(P_n⊗I)·C
+            # = (C^{-1}(P_1⊗I)C)·...·(C^{-1}(P_n⊗I)C)
+            # For the Pauli on the kth qubit P_k. The conjugation is calculated as following.
+            #   Puali X_k's conjugation is from the destabilzer table;
+            #   Puali Z_k's conjugation is from the stabilzer table;
+            #   Puali Y_k's conjugation is calcluated according to Y = iXZ. E.g., for the kth qubit,
+            #     C^{-1}·Y_k⊗I·C = C^{-1}·(iX_k⊗I·Z_k⊗I)·C = i (C^{-1}·X_k⊗I·C)·(C^{-1}·Z_k⊗I·C)
+            for qid, qubit in enumerate(op.qubits):
+                pauli = ps.get(qubit)
+                match pauli:
+                    case None:
+                        continue
+                    case pauli_gates.X:
+                        conjugated *= tableau.destabilizers()[qid]
+                    case pauli_gates.Z:
+                        conjugated *= tableau.stabilizers()[qid]
+                    case pauli_gates.Y:
+                        conjugated *= (
+                            1j
+                            * tableau.destabilizers()[qid]  # conj X first
+                            * tableau.stabilizers()[qid]  # then conj Z
+                        )
+            ps = remain * conjugated.on(*op.qubits)
+        return ps
 
     def after(self, ops: 'cirq.OP_TREE') -> 'cirq.PauliString':
         r"""Determines the equivalent pauli string after some operations.

cirq-core/cirq/ops/pauli_string_test.py

@@ -57,6 +57,22 @@ def _small_sample_qubit_pauli_maps():
     yield {qubits[0]: cirq.Z, qubits[1]: cirq.X, qubits[2]: cirq.Y}
 
 
+def assert_conjugation(
+    input_ps: cirq.PauliString, ops: cirq.OP_TREE, expected: cirq.PauliString | None
+):
+    conjugation = input_ps.conjugated_by(ops)
+    if expected is not None:
+        assert conjugation == expected
+    else:  # Compares the unitary of the conjugation result and the expected unitary.
+        op_list = list(cirq.flatten_to_ops(ops))
+        qubits_of_clifford = [q for op in op_list for q in op.qubits]
+        clifford = cirq.CliffordGate.from_op_list(op_list, qubits_of_clifford)
+        actual_unitary = cirq.unitary(conjugation.dense(qubits_of_clifford))
+        c = cirq.unitary(clifford)
+        expected_unitary = np.conj(c.T) @ cirq.unitary(input_ps.dense(qubits_of_clifford)) @ c
+        assert np.allclose(actual_unitary, expected_unitary, atol=1e-8)
+
+
 def test_eq_ne_hash():
     q0, q1, q2 = _make_qubits(3)
     eq = cirq.testing.EqualsTester()
@@ -1381,13 +1397,24 @@ def test_pauli_string_expectation_from_state_vector_mixed_state_linearity():
 def test_conjugated_by_normal_gates():
     a = cirq.LineQubit(0)
 
-    assert cirq.X(a).conjugated_by(cirq.H(a)) == cirq.Z(a)
-    assert cirq.Y(a).conjugated_by(cirq.H(a)) == -cirq.Y(a)
-    assert cirq.Z(a).conjugated_by(cirq.H(a)) == cirq.X(a)
+    assert_conjugation(cirq.X(a), cirq.H(a), cirq.Z(a))
+    assert_conjugation(cirq.Y(a), cirq.H(a), -cirq.Y(a))
+    assert_conjugation(cirq.Z(a), cirq.H(a), cirq.X(a))
+
+    assert_conjugation(cirq.X(a), cirq.S(a), -cirq.Y(a))
+    assert_conjugation(cirq.Y(a), cirq.S(a), cirq.X(a))
+    assert_conjugation(cirq.Z(a), cirq.S(a), cirq.Z(a))
+
+    clifford_op = cirq.PhasedXZGate(axis_phase_exponent=0.25, x_exponent=-1, z_exponent=0).on(a)
+    assert_conjugation(cirq.X(a), clifford_op, cirq.Y(a))
+    assert_conjugation(cirq.Y(a), clifford_op, cirq.X(a))
+    assert_conjugation(cirq.Z(a), clifford_op, -cirq.Z(a))
+
+
+def test_conjugated_by_op_gate_of_clifford_gate_type():
+    a = cirq.LineQubit(0)
 
-    assert cirq.X(a).conjugated_by(cirq.S(a)) == -cirq.Y(a)
-    assert cirq.Y(a).conjugated_by(cirq.S(a)) == cirq.X(a)
-    assert cirq.Z(a).conjugated_by(cirq.S(a)) == cirq.Z(a)
+    assert_conjugation(cirq.X(a), cirq.CliffordGate.from_op_list([cirq.H(a)], [a]).on(a), cirq.Z(a))
 
 
 def test_dense():
@@ -1430,16 +1457,25 @@ def test_conjugated_by_incorrectly_powered_cliffords():
         cirq.ZZ(a, b),
     ]
     for c in cliffords:
-        with pytest.raises(TypeError, match='not a known Clifford'):
+        with pytest.raises(
+            ValueError,
+            match='Clifford Gate can only be constructed from the operations'
+            ' that has stabilizer effect.',
+        ):
             _ = p.conjugated_by(c**0.1)
-        with pytest.raises(TypeError, match='not a known Clifford'):
+        with pytest.raises(
+            ValueError,
+            match='Clifford Gate can only be constructed from the operations'
+            ' that has stabilizer effect.',
+        ):
             _ = p.conjugated_by(c ** sympy.Symbol('t'))
 
 
 def test_conjugated_by_global_phase():
+    """Global phase gate preserves PauliString."""
     a = cirq.LineQubit(0)
-    assert cirq.X(a).conjugated_by(cirq.global_phase_operation(1j)) == cirq.X(a)
-    assert cirq.Z(a).conjugated_by(cirq.global_phase_operation(np.exp(1.1j))) == cirq.Z(a)
+    assert_conjugation(cirq.X(a), cirq.global_phase_operation(1j), cirq.X(a))
+    assert_conjugation(cirq.X(a), cirq.global_phase_operation(np.exp(1.1j)), cirq.X(a))
 
     class DecomposeGlobal(cirq.Gate):
         def num_qubits(self):
@@ -1448,7 +1484,7 @@ def num_qubits(self):
         def _decompose_(self, qubits):
             yield cirq.global_phase_operation(1j)
 
-    assert cirq.X(a).conjugated_by(DecomposeGlobal().on(a)) == cirq.X(a)
+    assert_conjugation(cirq.X(a), DecomposeGlobal().on(a), cirq.X(a))
 
 
 def test_conjugated_by_composite_with_disjoint_sub_gates():
@@ -1461,8 +1497,10 @@ def num_qubits(self):
         def _decompose_(self, qubits):
             yield cirq.H(qubits[1])
 
-    assert cirq.X(a).conjugated_by(DecomposeDisjoint().on(a, b)) == cirq.X(a)
-    assert cirq.X(a).pass_operations_over([DecomposeDisjoint().on(a, b)]) == cirq.X(a)
+    for g1 in [cirq.X, cirq.Y]:
+        for g2 in [cirq.X, cirq.Y]:
+            ps = g1(a) * g2(b)
+            assert ps.conjugated_by(DecomposeDisjoint().on(a, b)) == ps.conjugated_by(cirq.H(b))
 
 
 def test_conjugated_by_clifford_composite():
@@ -1477,16 +1515,16 @@ def _decompose_(self, qubits):
             yield cirq.SWAP(qubits[2], qubits[3])
 
     a, b, c, d = cirq.LineQubit.range(4)
-    p = cirq.X(a) * cirq.Z(b)
+    ps = cirq.X(a) * cirq.Z(b)
     u = UnknownGate()
-    assert p.conjugated_by(u(a, b, c, d)) == cirq.Z(a) * cirq.X(b)
+    assert_conjugation(ps, u(a, b, c, d), cirq.Z(a) * cirq.X(b))
 
 
 def test_conjugated_by_move_into_uninvolved():
     a, b, c, d = cirq.LineQubit.range(4)
-    p = cirq.X(a) * cirq.Z(b)
-    assert p.conjugated_by([cirq.SWAP(c, d), cirq.SWAP(b, c)]) == cirq.X(a) * cirq.Z(d)
-    assert p.conjugated_by([cirq.SWAP(b, c), cirq.SWAP(c, d)]) == cirq.X(a) * cirq.Z(c)
+    ps = cirq.X(a) * cirq.Z(b)
+    assert_conjugation(ps, [cirq.SWAP(c, d), cirq.SWAP(b, c)], cirq.X(a) * cirq.Z(d))
+    assert_conjugation(ps, [cirq.SWAP(b, c), cirq.SWAP(c, d)], cirq.X(a) * cirq.Z(c))
 
 
 def test_conjugated_by_common_single_qubit_gates():
@@ -1508,21 +1546,13 @@ def test_conjugated_by_common_single_qubit_gates():
     single_qubit_gates = [g**i for i in range(4) for g in base_single_qubit_gates]
     for p in [cirq.X, cirq.Y, cirq.Z]:
         for g in single_qubit_gates:
-            assert p.on(a).conjugated_by(g.on(b)) == p.on(a)
-
-            actual = cirq.unitary(p.on(a).conjugated_by(g.on(a)))
-            u = cirq.unitary(g)
-            expected = np.conj(u.T) @ cirq.unitary(p) @ u
-            assert cirq.allclose_up_to_global_phase(actual, expected, atol=1e-8)
+            # pauli gate on a, clifford on b: pauli gate preserves.
+            assert_conjugation(p(a), g(b), p(a))
+            # pauli gate on a, clifford on a: check conjugation in matrices.
+            assert_conjugation(p(a), g(a), None)
 
 
 def test_conjugated_by_common_two_qubit_gates():
-    class OrderSensitiveGate(cirq.Gate):
-        def num_qubits(self):
-            return 2
-
-        def _decompose_(self, qubits):
-            return [cirq.Y(qubits[0]) ** -0.5, cirq.CNOT(*qubits)]
 
     a, b, c, d = cirq.LineQubit.range(4)
     two_qubit_gates = [
@@ -1541,34 +1571,25 @@ def _decompose_(self, qubits):
         cirq.YY**-0.5,
         cirq.ZZ**-0.5,
     ]
-    two_qubit_gates.extend([OrderSensitiveGate()])
     for p1 in [cirq.I, cirq.X, cirq.Y, cirq.Z]:
         for p2 in [cirq.I, cirq.X, cirq.Y, cirq.Z]:
             pd = cirq.DensePauliString([p1, p2])
-            p = pd.sparse()
+            p = pd.sparse([a, b])
             for g in two_qubit_gates:
-                assert p.conjugated_by(g.on(c, d)) == p
-
-                actual = cirq.unitary(p.conjugated_by(g.on(a, b)).dense([a, b]))
-                u = cirq.unitary(g)
-                expected = np.conj(u.T) @ cirq.unitary(pd) @ u
-                np.testing.assert_allclose(actual, expected, atol=1e-8)
+                # pauli_string on (a,b), clifford on (c,d): pauli_string preserves.
+                assert_conjugation(p, g(c, d), p)
+                # pauli_string on (a,b), clifford on (a,b): compare unitaries of
+                # the conjugated_by and actual matrix conjugation.
+                assert_conjugation(p, g.on(a, b), None)
 
 
 def test_conjugated_by_ordering():
-    class OrderSensitiveGate(cirq.Gate):
-        def num_qubits(self):
-            return 2
-
-        def _decompose_(self, qubits):
-            return [cirq.Y(qubits[0]) ** -0.5, cirq.CNOT(*qubits)]
-
+    """Tests .conjugated_by([op1, op2]) == .conjugated_by(op2).conjugated_by(op1)"""
     a, b = cirq.LineQubit.range(2)
     inp = cirq.Z(b)
-    out1 = inp.conjugated_by(OrderSensitiveGate().on(a, b))
-    out2 = inp.conjugated_by([cirq.H(a), cirq.CNOT(a, b)])
-    out3 = inp.conjugated_by(cirq.CNOT(a, b)).conjugated_by(cirq.H(a))
-    assert out1 == out2 == out3 == cirq.X(a) * cirq.Z(b)
+    out1 = inp.conjugated_by([cirq.H(a), cirq.CNOT(a, b)])
+    out2 = inp.conjugated_by(cirq.CNOT(a, b)).conjugated_by(cirq.H(a))
+    assert out1 == out2 == cirq.X(a) * cirq.Z(b)
 
 
 def test_pass_operations_over_ordering():

cirq-core/cirq/ops/swap_gates.py

@@ -219,6 +219,11 @@ def _eigen_components(self) -> List[Tuple[float, np.ndarray]]:
         ]
         # yapf: enable
 
+    def _has_stabilizer_effect_(self) -> Optional[bool]:
+        if self._is_parameterized_():
+            return None
+        return self.exponent % 1 == 0
+
     def _decompose_(self, qubits):
         a, b = qubits