diff --git a/src/qibo/gates/abstract.py b/src/qibo/gates/abstract.py
index 2cc99fa4b5..6476883d69 100644
--- a/src/qibo/gates/abstract.py
+++ b/src/qibo/gates/abstract.py
@@ -48,9 +48,9 @@ def __init__(self):
 
         self.clifford = False
         self.unitary = False
-        self._target_qubits = tuple()
-        self._control_qubits = set()
-        self._parameters = tuple()
+        self._target_qubits = ()
+        self._control_qubits = ()
+        self._parameters = ()
         config.ALLOW_SWITCHERS = False
 
         self.symbolic_parameters = {}
@@ -158,13 +158,13 @@ def _set_target_qubits(self, qubits: Sequence[int]):
 
     def _set_control_qubits(self, qubits: Sequence[int]):
         """Helper method for setting control qubits."""
-        self._control_qubits = set(qubits)
-        if len(self._control_qubits) != len(qubits):
+        if len(set(qubits)) != len(qubits):
             repeated = self._find_repeated(qubits)
             raise_error(
                 ValueError,
                 f"Control qubit {repeated} was given twice for gate {self.__class__.__name__}.",
             )
+        self._control_qubits = qubits
 
     @target_qubits.setter
     def target_qubits(self, qubits: Sequence[int]):
@@ -204,11 +204,12 @@ def _find_repeated(qubits: Sequence[int]) -> int:
     def _check_control_target_overlap(self):
         """Checks that there are no qubits that are both target and
         controls."""
-        common = set(self._target_qubits) & self._control_qubits
+        control_and_target = self._control_qubits + self._target_qubits
+        common = len(set(control_and_target)) != len(control_and_target)
         if common:
             raise_error(
                 ValueError,
-                f"{common} qubits are both targets and controls "
+                f"{set(self._target_qubits) & set(self._control_qubits)} qubits are both targets and controls "
                 + f"for gate {self.__class__.__name__}.",
             )
 
diff --git a/src/qibo/quantum_info/random_ensembles.py b/src/qibo/quantum_info/random_ensembles.py
index e14db95f35..26e6932f1b 100644
--- a/src/qibo/quantum_info/random_ensembles.py
+++ b/src/qibo/quantum_info/random_ensembles.py
@@ -2,6 +2,7 @@
 
 import math
 import warnings
+from functools import cache
 from typing import Optional, Union
 
 import numpy as np
@@ -1101,6 +1102,21 @@ def _sample_from_quantum_mallows_distribution(nqubits: int, local_state):
     return hadamards, permutations
 
 
+@cache
+def _create_S(q):
+    return gates.S(q)
+
+
+@cache
+def _create_CZ(cq, tq):
+    return gates.CZ(cq, tq)
+
+
+@cache
+def _create_CNOT(cq, tq):
+    return gates.CNOT(cq, tq)
+
+
 def _operator_from_hadamard_free_group(
     gamma_matrix, delta_matrix, density_matrix: bool = False, pauli_operator=None
 ):
@@ -1138,19 +1154,19 @@ def _operator_from_hadamard_free_group(
     if pauli_operator is not None:
         circuit += pauli_operator
 
-    for qubit, gamma in enumerate(np.diag(gamma_matrix)):
-        if gamma == 1:
-            circuit.add(gates.S(qubit))
+    idx = np.tril_indices(nqubits, k=-1)
+    gamma_ones = gamma_matrix[idx].nonzero()[0]
+    delta_ones = delta_matrix[idx].nonzero()[0]
 
-    for j in range(nqubits):
-        for k in range(j + 1, nqubits):
-            if gamma_matrix[k, j] == 1:
-                circuit.add(gates.CZ(j, k))
+    S_gates = [_create_S(q) for q in np.diag(gamma_matrix).nonzero()[0]]
+    CZ_gates = [
+        _create_CZ(cq, tq) for cq, tq in zip(idx[1][gamma_ones], idx[0][gamma_ones])
+    ]
+    CNOT_gates = [
+        _create_CNOT(cq, tq) for cq, tq in zip(idx[1][delta_ones], idx[0][delta_ones])
+    ]
 
-    for j in range(nqubits):
-        for k in range(j + 1, nqubits):
-            if delta_matrix[k, j] == 1:
-                circuit.add(gates.CNOT(j, k))
+    circuit.add(S_gates + CZ_gates + CNOT_gates)
 
     return circuit