Extend Clifford protocol for 2+ qubit operations through analytical c…

…heck and falling back to decompose (#6332)

cirq-core/cirq/protocols/has_stabilizer_effect_protocol.py

@@ -15,8 +15,18 @@
 from typing import Any, Optional
 
 from cirq.ops.clifford_gate import SingleQubitCliffordGate
+from cirq.ops.dense_pauli_string import DensePauliString
+from cirq._import import LazyLoader
+import cirq.protocols.unitary_protocol as unitary_protocol
+import cirq.protocols.has_unitary_protocol as has_unitary_protocol
+import cirq.protocols.qid_shape_protocol as qid_shape_protocol
+import cirq.protocols.decompose_protocol as decompose_protocol
 
-from cirq import protocols
+pauli_string_decomposition = LazyLoader(
+    "pauli_string_decomposition",
+    globals(),
+    "cirq.transformers.analytical_decompositions.pauli_string_decomposition",
+)
 
 
 def has_stabilizer_effect(val: Any) -> bool:
@@ -29,6 +39,7 @@ def has_stabilizer_effect(val: Any) -> bool:
         _strat_has_stabilizer_effect_from_has_stabilizer_effect,
         _strat_has_stabilizer_effect_from_gate,
         _strat_has_stabilizer_effect_from_unitary,
+        _strat_has_stabilizer_effect_from_decompose,
     ]
     for strat in strats:
         result = strat(val)
@@ -62,9 +73,44 @@ def _strat_has_stabilizer_effect_from_unitary(val: Any) -> Optional[bool]:
     2x2 unitaries.
     """
     # Do not try this strategy if there is no unitary or if the number of
-    # qubits is not 1 since that would be expensive.
-    qid_shape = protocols.qid_shape(val, default=None)
-    if qid_shape is None or len(qid_shape) != 1 or not protocols.has_unitary(val):
+    # qubits is greater than 3 since that would be expensive.
+    qid_shape = qid_shape_protocol.qid_shape(val, default=None)
+    if (
+        qid_shape is None
+        or len(qid_shape) > 3
+        or qid_shape != (2,) * len(qid_shape)
+        or not has_unitary_protocol.has_unitary(val)
+    ):
         return None
-    unitary = protocols.unitary(val)
-    return SingleQubitCliffordGate.from_unitary(unitary) is not None
+    unitary = unitary_protocol.unitary(val)
+    if len(qid_shape) == 1:
+        return SingleQubitCliffordGate.from_unitary(unitary) is not None
+
+    # Check if the action of the unitary on each single qubit pauli string leads to a pauli product.
+    # Source: https://quantumcomputing.stackexchange.com/a/13158
+    for q_idx in range(len(qid_shape)):
+        for g in 'XZ':
+            pauli_string = ['I'] * len(qid_shape)
+            pauli_string[q_idx] = g
+            ps = DensePauliString(pauli_string)
+            p = ps._unitary_()
+            if not pauli_string_decomposition.unitary_to_pauli_string(
+                (unitary @ p @ unitary.T.conj())
+            ):
+                return False
+    return True
+
+
+def _strat_has_stabilizer_effect_from_decompose(val: Any) -> Optional[bool]:
+    qid_shape = qid_shape_protocol.qid_shape(val, default=None)
+    if qid_shape is None or len(qid_shape) <= 3:
+        return None
+
+    decomposition = decompose_protocol.decompose_once(val, default=None)
+    if decomposition is None:
+        return None
+    for op in decomposition:
+        res = has_stabilizer_effect(op)
+        if not res:
+            return res
+    return True

cirq-core/cirq/protocols/has_stabilizer_effect_protocol_test.py

@@ -92,6 +92,10 @@ def __init__(self, unitary):
     def _unitary_(self):
         return self.unitary
 
+    @property
+    def qubits(self):
+        return cirq.LineQubit.range(self.unitary.shape[0].bit_length() - 1)
+
 
 def test_inconclusive():
     assert not cirq.has_stabilizer_effect(object())
@@ -125,9 +129,20 @@ def test_via_unitary():
     op3 = OpWithUnitary(np.array([[1, 0], [0, np.sqrt(1j)]]))
     assert not cirq.has_stabilizer_effect(op3)
 
+    # 2+ qubit cliffords
+    assert cirq.has_stabilizer_effect(cirq.CNOT)
+    assert cirq.has_stabilizer_effect(cirq.XX)
+    assert cirq.has_stabilizer_effect(cirq.ZZ)
+
+    # Non Cliffords
+    assert not cirq.has_stabilizer_effect(cirq.T)
+    assert not cirq.has_stabilizer_effect(cirq.CCNOT)
+    assert not cirq.has_stabilizer_effect(cirq.CCZ)
+
 
-def test_via_unitary_not_supported():
-    # Unitaries larger than 2x2 are not yet supported.
-    op = OpWithUnitary(cirq.unitary(cirq.CNOT))
-    assert not cirq.has_stabilizer_effect(op)
-    assert not cirq.has_stabilizer_effect(op)
+def test_via_decompose():
+    assert cirq.has_stabilizer_effect(cirq.Circuit(cirq.H.on_each(cirq.LineQubit.range(4))))
+    assert not cirq.has_stabilizer_effect(cirq.Circuit(cirq.T.on_each(cirq.LineQubit.range(4))))
+    assert not cirq.has_stabilizer_effect(
+        OpWithUnitary(cirq.unitary(cirq.Circuit(cirq.T.on_each(cirq.LineQubit.range(4)))))
+    )

cirq-ft/cirq_ft/infra/t_complexity_protocol.py

@@ -81,7 +81,7 @@ def _is_clifford_or_t(stc: Any, fail_quietly: bool) -> Optional[TComplexity]:
     if isinstance(stc, cirq.ClassicallyControlledOperation):
         stc = stc.without_classical_controls()
 
-    if cirq.has_stabilizer_effect(stc):
+    if cirq.num_qubits(stc) <= 2 and cirq.has_stabilizer_effect(stc):
         # Clifford operation.
         return TComplexity(clifford=1)