make SingleQubitCliffordGate immutable singletons and use it in qubit_characterizations for a 37% speedup (#6392)

SingleQubitCliffordGates represents the 24 gates in the single qubit clifford group. Some of the operations this class implements are expensive computation but at the same time are properties of each of the 24 operators.

turning those expensive computations into cached properties is benificial for performance but for that the objects need to be immutable.

one of the places that heavily relies on single qubit cliffords is `qubit_characterizations.py` which implemented the single qubit clifford algebra from scratch by falling on to the matrix representation and doing matrix multiplication and inversion (for computing the adjoint) this led to a bottleneck while creating circuits (e.g. for example _create_parallel_rb_circuit for 50 qubits and a 1000 gates takes $3.886s$). using SingleQubitCliffordGates instead and using `merged_with` operation which maps two cliffords onto the clifford equaivalent to their composition leads to $2.148s$, with most of those 2s spent in `Moment.__init__` which will be the target of the next PR.

I also made the 24 cliffords singleton, since there is no point in creating new object which won't have any of the cached properties.
 
part of #6389

Author: NoureldinYosri

cirq-core/cirq/experiments/qubit_characterizations.py

@@ -14,6 +14,7 @@
 
 import dataclasses
 import itertools
+import functools
 
 from typing import (
     Any,
@@ -58,11 +59,11 @@ class Cliffords:
         s1_y
     """
 
-    c1_in_xy: List[List[ops.Gate]]
-    c1_in_xz: List[List[ops.Gate]]
-    s1: List[List[ops.Gate]]
-    s1_x: List[List[ops.Gate]]
-    s1_y: List[List[ops.Gate]]
+    c1_in_xy: List[List[ops.SingleQubitCliffordGate]]
+    c1_in_xz: List[List[ops.SingleQubitCliffordGate]]
+    s1: List[List[ops.SingleQubitCliffordGate]]
+    s1_x: List[List[ops.SingleQubitCliffordGate]]
+    s1_y: List[List[ops.SingleQubitCliffordGate]]
 
 
 class RandomizedBenchMarkResult:
@@ -299,17 +300,14 @@ def parallel_single_qubit_randomized_benchmarking(
         A dictionary from qubits to RandomizedBenchMarkResult objects.
     """
 
-    cliffords = _single_qubit_cliffords()
-    c1 = cliffords.c1_in_xy if use_xy_basis else cliffords.c1_in_xz
-    clifford_mats = np.array([_gate_seq_to_mats(gates) for gates in c1])
+    clifford_group = _single_qubit_cliffords()
+    c1 = clifford_group.c1_in_xy if use_xy_basis else clifford_group.c1_in_xz
 
     # create circuits
     circuits_all: List['cirq.AbstractCircuit'] = []
     for num_cliffords in num_clifford_range:
         for _ in range(num_circuits):
-            circuits_all.append(
-                _create_parallel_rb_circuit(qubits, num_cliffords, c1, clifford_mats)
-            )
+            circuits_all.append(_create_parallel_rb_circuit(qubits, num_cliffords, c1))
 
     # run circuits
     results = sampler.run_batch(circuits_all, repetitions=repetitions)
@@ -562,11 +560,9 @@ def _measurement(two_qubit_circuit: circuits.Circuit) -> np.ndarray:
 
 
 def _create_parallel_rb_circuit(
-    qubits: Iterator['cirq.Qid'], num_cliffords: int, c1: list, clifford_mats: np.ndarray
+    qubits: Iterator['cirq.Qid'], num_cliffords: int, c1: list
 ) -> 'cirq.Circuit':
-    circuits_to_zip = [
-        _random_single_q_clifford(qubit, num_cliffords, c1, clifford_mats) for qubit in qubits
-    ]
+    circuits_to_zip = [_random_single_q_clifford(qubit, num_cliffords, c1) for qubit in qubits]
     circuit = circuits.Circuit.zip(*circuits_to_zip)
     return circuits.Circuit.from_moments(*circuit, ops.measure_each(*qubits))
 
@@ -612,16 +608,12 @@ def _two_qubit_clifford_matrices(
 
 
 def _random_single_q_clifford(
-    qubit: 'cirq.Qid',
-    num_cfds: int,
-    cfds: Sequence[Sequence['cirq.Gate']],
-    cfd_matrices: np.ndarray,
+    qubit: 'cirq.Qid', num_cfds: int, cfds: Sequence[Sequence['cirq.Gate']]
 ) -> 'cirq.Circuit':
     clifford_group_size = 24
     gate_ids = list(np.random.choice(clifford_group_size, num_cfds))
     gate_sequence = [gate for gate_id in gate_ids for gate in cfds[gate_id]]
-    idx = _find_inv_matrix(_gate_seq_to_mats(gate_sequence), cfd_matrices)
-    gate_sequence.extend(cfds[idx])
+    gate_sequence.append(_reduce_gate_seq(gate_sequence) ** -1)
     circuit = circuits.Circuit(gate(qubit) for gate in gate_sequence)
     return circuit
 
@@ -681,11 +673,13 @@ def _matrix_bar_plot(
         ax.set_title(title)
 
 
-def _gate_seq_to_mats(gate_seq: Sequence['cirq.Gate']) -> np.ndarray:
-    mat_rep = protocols.unitary(gate_seq[0])
+def _reduce_gate_seq(
+    gate_seq: Sequence[ops.SingleQubitCliffordGate],
+) -> ops.SingleQubitCliffordGate:
+    cur = gate_seq[0]
     for gate in gate_seq[1:]:
-        mat_rep = np.dot(protocols.unitary(gate), mat_rep)
-    return mat_rep
+        cur = cur.merged_with(gate)
+    return cur
 
 
 def _two_qubit_clifford(
@@ -793,11 +787,16 @@ def _single_qubit_gates(
         yield gate(qubit)
 
 
+@functools.cache
 def _single_qubit_cliffords() -> Cliffords:
-    X, Y, Z = ops.X, ops.Y, ops.Z
+    X, Y, Z = (
+        ops.SingleQubitCliffordGate.X,
+        ops.SingleQubitCliffordGate.Y,
+        ops.SingleQubitCliffordGate.Z,
+    )
 
-    c1_in_xy: List[List['cirq.Gate']] = []
-    c1_in_xz: List[List['cirq.Gate']] = []
+    c1_in_xy: List[List[ops.SingleQubitCliffordGate]] = []
+    c1_in_xz: List[List[ops.SingleQubitCliffordGate]] = []
 
     for phi_0, phi_1 in itertools.product([1.0, 0.5, -0.5], [0.0, 0.5, -0.5]):
         c1_in_xy.append([X**phi_0, Y**phi_1])
@@ -820,8 +819,20 @@ def _single_qubit_cliffords() -> Cliffords:
     for z0, x, z1 in phi_xz:
         c1_in_xz.append([Z**z0, X**x, Z**z1])
 
-    s1: List[List['cirq.Gate']] = [[X**0.0], [Y**0.5, X**0.5], [X**-0.5, Y**-0.5]]
-    s1_x: List[List['cirq.Gate']] = [[X**0.5], [X**0.5, Y**0.5, X**0.5], [Y**-0.5]]
-    s1_y: List[List['cirq.Gate']] = [[Y**0.5], [X**-0.5, Y**-0.5, X**0.5], [Y, X**0.5]]
+    s1: List[List[ops.SingleQubitCliffordGate]] = [
+        [X**0.0],
+        [Y**0.5, X**0.5],
+        [X**-0.5, Y**-0.5],
+    ]
+    s1_x: List[List[ops.SingleQubitCliffordGate]] = [
+        [X**0.5],
+        [X**0.5, Y**0.5, X**0.5],
+        [Y**-0.5],
+    ]
+    s1_y: List[List[ops.SingleQubitCliffordGate]] = [
+        [Y**0.5],
+        [X**-0.5, Y**-0.5, X**0.5],
+        [Y, X**0.5],
+    ]
 
     return Cliffords(c1_in_xy, c1_in_xz, s1, s1_x, s1_y)

cirq-core/cirq/experiments/qubit_characterizations_test.py

@@ -75,9 +75,32 @@ def check_distinct(unitaries):
 
     # Check that XZ decomposition has at most one X gate per clifford.
     for gates in cliffords.c1_in_xz:
-        num_x = len([gate for gate in gates if isinstance(gate, cirq.XPowGate)])
-        num_z = len([gate for gate in gates if isinstance(gate, cirq.ZPowGate)])
-        assert num_x + num_z == len(gates)
+        num_i = len([gate for gate in gates if gate == cirq.ops.SingleQubitCliffordGate.I])
+        num_x = len(
+            [
+                gate
+                for gate in gates
+                if gate
+                in (
+                    cirq.ops.SingleQubitCliffordGate.X,
+                    cirq.ops.SingleQubitCliffordGate.X_sqrt,
+                    cirq.ops.SingleQubitCliffordGate.X_nsqrt,
+                )
+            ]
+        )
+        num_z = len(
+            [
+                gate
+                for gate in gates
+                if gate
+                in (
+                    cirq.ops.SingleQubitCliffordGate.Z,
+                    cirq.ops.SingleQubitCliffordGate.Z_sqrt,
+                    cirq.ops.SingleQubitCliffordGate.Z_nsqrt,
+                )
+            ]
+        )
+        assert num_x + num_z + num_i == len(gates)
         assert num_x <= 1
 
 

cirq-core/cirq/ops/clifford_gate.py

@@ -14,11 +14,13 @@
 
 from typing import Any, Dict, List, Optional, Sequence, Tuple, TYPE_CHECKING, Union
 
-
+import functools
+from dataclasses import dataclass
 import numpy as np
 
 from cirq import protocols, value, linalg, qis
 from cirq._import import LazyLoader
+from cirq._compat import cached_property, cached_method
 from cirq.ops import common_gates, named_qubit, raw_types, pauli_gates, phased_x_z_gate
 from cirq.ops.pauli_gates import Pauli
 from cirq.type_workarounds import NotImplementedType
@@ -356,6 +358,8 @@ def _get_sqrt_map(
 class CliffordGate(raw_types.Gate, CommonCliffordGates):
     """Clifford rotation for N-qubit."""
 
+    _clifford_tableau: qis.CliffordTableau
+
     def __init__(self, *, _clifford_tableau: qis.CliffordTableau) -> None:
         # We use the Clifford tableau to represent a Clifford gate.
         # It is crucial to note that the meaning of tableau here is different
@@ -376,7 +380,7 @@ def __init__(self, *, _clifford_tableau: qis.CliffordTableau) -> None:
         # more precisely the conjugate transformation of ZI by this gate, becomes -ZI.
         # (Note the real clifford tableau has to satify the Symplectic property.
         # here is just for illustration)
-        self._clifford_tableau = _clifford_tableau.copy()
+        object.__setattr__(self, '_clifford_tableau', _clifford_tableau.copy())
 
     @property
     def clifford_tableau(self):
@@ -399,6 +403,12 @@ def _has_stabilizer_effect_(self) -> Optional[bool]:
     def __pow__(self, exponent) -> 'CliffordGate':
         if exponent == -1:
             return CliffordGate.from_clifford_tableau(self.clifford_tableau.inverse())
+        if exponent == 0:
+            return CliffordGate.from_clifford_tableau(
+                qis.CliffordTableau(num_qubits=self._num_qubits_())
+            )
+        if exponent == 1:
+            return self
         if exponent > 0 and int(exponent) == exponent:
             base_tableau = self.clifford_tableau.copy()
             for _ in range(int(exponent) - 1):
@@ -457,6 +467,7 @@ def _act_on_(
         return NotImplemented
 
 
+@dataclass(frozen=True, init=False, eq=False, repr=False)
 @value.value_equality(manual_cls=True)
 class SingleQubitCliffordGate(CliffordGate):
     """Any single qubit Clifford rotation."""
@@ -468,6 +479,7 @@ def _num_qubits_(self):
         return 1
 
     @staticmethod
+    @functools.cache
     def from_clifford_tableau(tableau: qis.CliffordTableau) -> 'SingleQubitCliffordGate':
         if not isinstance(tableau, qis.CliffordTableau):
             raise ValueError('Input argument has to be a CliffordTableau instance.')
@@ -679,6 +691,10 @@ def to_phased_xz_gate(self) -> phased_x_z_gate.PhasedXZGate:
             * {middle point of xyz in 4 Quadrant} * 120 is [[0, 1], [1, 1]]
             * {middle point of xyz in 4 Quadrant} * 240 is [[1, 1], [1, 0]]
         """
+        return self._to_phased_xz_gate
+
+    @cached_property
+    def _to_phased_xz_gate(self) -> phased_x_z_gate.PhasedXZGate:
         x_to_flip, z_to_flip = self.clifford_tableau.rs
         flip_index = int(z_to_flip) * 2 + x_to_flip
         a, x, z = 0.0, 0.0, 0.0
@@ -716,7 +732,7 @@ def to_phased_xz_gate(self) -> phased_x_z_gate.PhasedXZGate:
             z = -0.5 if x_to_flip else 0.5
         return phased_x_z_gate.PhasedXZGate(x_exponent=x, z_exponent=z, axis_phase_exponent=a)
 
-    def __pow__(self, exponent) -> 'SingleQubitCliffordGate':
+    def __pow__(self, exponent: Union[float, int]) -> 'SingleQubitCliffordGate':
         # First to check if we can get the sqrt and negative sqrt Clifford.
         if self._get_sqrt_map().get(exponent, None):
             pow_gate = self._get_sqrt_map()[exponent].get(self, None)
@@ -761,6 +777,7 @@ def commutes_with_pauli(self, pauli: Pauli) -> bool:
         to, flip = self.pauli_tuple(pauli)
         return to == pauli and not flip
 
+    @cached_method
     def merged_with(self, second: 'SingleQubitCliffordGate') -> 'SingleQubitCliffordGate':
         """Returns a SingleQubitCliffordGate such that the circuits
             --output-- and --self--second--
@@ -773,6 +790,10 @@ def _has_unitary_(self) -> bool:
         return True
 
     def _unitary_(self) -> np.ndarray:
+        return self._unitary
+
+    @cached_property
+    def _unitary(self) -> np.ndarray:
         mat = np.eye(2)
         qubit = named_qubit.NamedQubit('arbitrary')
         for op in protocols.decompose_once_with_qubits(self, (qubit,)):
@@ -787,6 +808,10 @@ def decompose_gate(self) -> Sequence['cirq.Gate']:
             clifford gate if applied in order. This decomposition agrees with
             cirq.unitary(self), including global phase.
         """
+        return self._decompose_gate
+
+    @cached_property
+    def _decompose_gate(self) -> Sequence['cirq.Gate']:
         if self == SingleQubitCliffordGate.H:
             return [common_gates.H]
         rotations = self.decompose_rotation()
@@ -802,6 +827,10 @@ def decompose_rotation(self) -> Sequence[Tuple[Pauli, int]]:
         Note that the combined unitary effect of these rotations may
         differ from cirq.unitary(self) by a global phase.
         """
+        return self._decompose_rotation
+
+    @cached_property
+    def _decompose_rotation(self) -> Sequence[Tuple[Pauli, int]]:
         x_rot = self.pauli_tuple(pauli_gates.X)
         y_rot = self.pauli_tuple(pauli_gates.Y)
         z_rot = self.pauli_tuple(pauli_gates.Z)
@@ -895,6 +924,10 @@ def _circuit_diagram_info_(
         )
 
     def _value_equality_values_(self):
+        return self._value_equality_values
+
+    @cached_property
+    def _value_equality_values(self):
         return self._clifford_tableau.matrix().tobytes() + self._clifford_tableau.rs.tobytes()
 
     def _value_equality_values_cls_(self):

cirq-core/cirq/qis/clifford_tableau.py

@@ -17,7 +17,7 @@
 import numpy as np
 
 from cirq import protocols
-from cirq._compat import proper_repr
+from cirq._compat import proper_repr, cached_method
 from cirq.qis import quantum_state_representation
 from cirq.value import big_endian_int_to_digits, linear_dict, random_state
 
@@ -652,3 +652,7 @@ def measure(
         self, axes: Sequence[int], seed: 'cirq.RANDOM_STATE_OR_SEED_LIKE' = None
     ) -> List[int]:
         return [self._measure(axis, random_state.parse_random_state(seed)) for axis in axes]
+
+    @cached_method
+    def __hash__(self) -> int:
+        return hash(self.matrix().tobytes() + self.rs.tobytes())