diff --git a/doc/source/code-examples/advancedexamples.rst b/doc/source/code-examples/advancedexamples.rst
index e6740b5141..572342f79f 100644
--- a/doc/source/code-examples/advancedexamples.rst
+++ b/doc/source/code-examples/advancedexamples.rst
@@ -289,12 +289,13 @@ The following gates support parameter setting:
 * :class:`qibo.gates.fSim`: Accepts a tuple of two parameters ``(theta, phi)``.
 * :class:`qibo.gates.GeneralizedfSim`: Accepts a tuple of two parameters
   ``(unitary, phi)``. Here ``unitary`` should be a unitary matrix given as an
-  array or ``tf.Tensor`` of shape ``(2, 2)``.
+  array or ``tf.Tensor`` of shape ``(2, 2)``. A ``torch.Tensor`` is required when using the pytorch backend.
 * :class:`qibo.gates.Unitary`: Accepts a single ``unitary`` parameter. This
-  should be an array or ``tf.Tensor`` of shape ``(2, 2)``.
+  should be an array or ``tf.Tensor`` of shape ``(2, 2)``. A ``torch.Tensor`` is required when using the pytorch backend.
 
 Note that a ``np.ndarray`` or a ``tf.Tensor`` may also be used in the place of
-a flat list. Using :meth:`qibo.models.circuit.Circuit.set_parameters` is more
+a flat list (``torch.Tensor`` is required when using the pytorch backend).
+Using :meth:`qibo.models.circuit.Circuit.set_parameters` is more
 efficient than recreating a new circuit with new parameter values. The inverse
 method :meth:`qibo.models.circuit.Circuit.get_parameters` is also available
 and returns a list, dictionary or flat list with the current parameter values
@@ -551,9 +552,9 @@ Here is a simple example using the Heisenberg XXZ model Hamiltonian:
 For more information on the available options of the ``vqe.minimize`` call we
 refer to the :ref:`Optimizers <Optimizers>` section of the documentation.
 Note that if the Stochastic Gradient Descent optimizer is used then the user
-has to use a backend based on tensorflow primitives and not the default custom
+has to use a backend based on tensorflow or pytorch primitives and not the default custom
 backend, as custom operators currently do not support automatic differentiation.
-To switch the backend one can do ``qibo.set_backend("tensorflow")``.
+To switch the backend one can do ``qibo.set_backend("tensorflow")`` or ``qibo.set_backend("pytorch")``.
 Check the :ref:`How to use automatic differentiation? <autodiff-example>`
 section for more details.
 
@@ -695,12 +696,13 @@ the model. For example the previous example would have to be modified as:
 How to use automatic differentiation?
 -------------------------------------
 
+The parameters of variational circuits can be optimized using the frameworks of
+Tensorflow or Pytorch.
+
 As a deep learning framework, Tensorflow supports
 `automatic differentiation <https://www.tensorflow.org/tutorials/customization/autodiff>`_.
-This can be used to optimize the parameters of variational circuits. For example
-the following script optimizes the parameters of two rotations so that the circuit
-output matches a target state using the fidelity as the corresponding loss
-function.
+The following script optimizes the parameters of two rotations so that the
+circuit output matches a target state using the fidelity as the corresponding loss function.
 
 Note that, as in the following example, the rotation angles have to assume real values
 to ensure the rotational gates are representing unitary operators.
@@ -777,6 +779,40 @@ that is supported by Tensorflow, such as defining
 and using the `Sequential model API <https://www.tensorflow.org/api_docs/python/tf/keras/Sequential>`_
 to train them.
 
+Similarly, Pytorch supports `automatic differentiation <https://pytorch.org/tutorials/beginner/basics/autogradqs_tutorFor%20example%20tial.html>`_.
+The following script optimizes the parameters of the variational circuit of the first example using the Pytorch framework.
+
+.. code-block:: python
+
+    import qibo
+    qibo.set_backend("pytorch")
+    import torch
+    from qibo import gates, models
+
+    # Optimization parameters
+    nepochs = 1000
+    optimizer = torch.optim.Adam
+    target_state = torch.ones(4, dtype=torch.complex128) / 2.0
+
+    # Define circuit ansatz
+    params = torch.tensor(
+        torch.rand(2, dtype=torch.float64), requires_grad=True
+    )
+    c = models.Circuit(2)
+    c.add(gates.RX(0, params[0]))
+    c.add(gates.RY(1, params[1]))
+
+    optimizer = optimizer([params])
+
+    for _ in range(nepochs):
+        optimizer.zero_grad()
+        c.set_parameters(params)
+        final_state = c().state()
+        fidelity = torch.abs(torch.sum(torch.conj(target_state) * final_state))
+        loss = 1 - fidelity
+        loss.backward()
+        optimizer.step()
+
 
 .. _noisy-example:
 
diff --git a/src/qibo/backends/npmatrices.py b/src/qibo/backends/npmatrices.py
index 69f0920722..723578bd0b 100644
--- a/src/qibo/backends/npmatrices.py
+++ b/src/qibo/backends/npmatrices.py
@@ -17,10 +17,6 @@ def __init__(self, dtype):
     def _cast(self, x, dtype):
         return self.np.array(x, dtype=dtype)
 
-    # This method is used to cast the parameters of the gates to the right type for other backends
-    def _cast_parameter(self, x):
-        return x
-
     @cached_property
     def H(self):
         return self._cast([[1, 1], [1, -1]], dtype=self.dtype) / math.sqrt(2)
@@ -66,34 +62,29 @@ def TDG(self):
         )
 
     def I(self, n=2):
-        return self._cast(self.np.eye(n), dtype=self.dtype)
+        return self.np.eye(n, dtype=self.dtype)
 
     def Align(self, delay, n=2):
-        return self._cast(self.I(n), dtype=self.dtype)
+        return self.I(n)
 
     def M(self):  # pragma: no cover
         raise_error(NotImplementedError)
 
     def RX(self, theta):
-        theta = self._cast_parameter(theta)
         cos = self.np.cos(theta / 2.0) + 0j
         isin = -1j * self.np.sin(theta / 2.0)
         return self._cast([[cos, isin], [isin, cos]], dtype=self.dtype)
 
     def RY(self, theta):
-        theta = self._cast_parameter(theta)
         cos = self.np.cos(theta / 2.0) + 0j
         sin = self.np.sin(theta / 2.0) + 0j
         return self._cast([[cos, -sin], [sin, cos]], dtype=self.dtype)
 
     def RZ(self, theta):
-        theta = self._cast_parameter(theta)
         phase = self.np.exp(0.5j * theta)
         return self._cast([[self.np.conj(phase), 0], [0, phase]], dtype=self.dtype)
 
     def PRX(self, theta, phi):
-        theta = self._cast_parameter(theta)
-        phi = self._cast_parameter(phi)
         cos = self.np.cos(theta / 2)
         sin = self.np.sin(theta / 2)
         exponent1 = -1.0j * self.np.exp(-1.0j * phi)
@@ -104,25 +95,20 @@ def PRX(self, theta, phi):
         )
 
     def GPI(self, phi):
-        phi = self._cast_parameter(phi)
         phase = self.np.exp(1.0j * phi)
         return self._cast([[0, self.np.conj(phase)], [phase, 0]], dtype=self.dtype)
 
     def GPI2(self, phi):
-        phi = self._cast_parameter(phi)
         phase = self.np.exp(1.0j * phi)
         return self._cast(
             [[1, -1.0j * self.np.conj(phase)], [-1.0j * phase, 1]], dtype=self.dtype
         ) / math.sqrt(2)
 
     def U1(self, theta):
-        theta = self._cast_parameter(theta)
         phase = self.np.exp(1j * theta)
         return self._cast([[1, 0], [0, phase]], dtype=self.dtype)
 
     def U2(self, phi, lam):
-        phi = self._cast_parameter(phi)
-        lam = self._cast_parameter(lam)
         eplus = self.np.exp(1j * (phi + lam) / 2.0)
         eminus = self.np.exp(1j * (phi - lam) / 2.0)
         return self._cast(
@@ -131,9 +117,6 @@ def U2(self, phi, lam):
         ) / math.sqrt(2)
 
     def U3(self, theta, phi, lam):
-        theta = self._cast_parameter(theta)
-        phi = self._cast_parameter(phi)
-        lam = self._cast_parameter(lam)
         cost = self.np.cos(theta / 2)
         sint = self.np.sin(theta / 2)
         eplus = self.np.exp(1j * (phi + lam) / 2.0)
@@ -147,8 +130,6 @@ def U3(self, theta, phi, lam):
         )
 
     def U1q(self, theta, phi):
-        theta = self._cast_parameter(theta)
-        phi = self._cast_parameter(phi)
         return self._cast(
             self.U3(theta, phi - math.pi / 2, math.pi / 2 - phi), dtype=self.dtype
         )
@@ -179,12 +160,12 @@ def CZ(self):
 
     @cached_property
     def CSX(self):
-        a = self._cast_parameter((1 + 1j) / 2)
-        b = self.np.conj(a)
+        a = (1 + 1j) / 2
+        b = (1 - 1j) / 2
         return self._cast(
             [
-                [1, 0, 0, 0],
-                [0, 1, 0, 0],
+                [1 + 0j, 0, 0, 0],
+                [0, 1 + 0j, 0, 0],
                 [0, 0, a, b],
                 [0, 0, b, a],
             ],
@@ -193,20 +174,19 @@ def CSX(self):
 
     @cached_property
     def CSXDG(self):
-        a = self._cast_parameter((1 - 1j) / 2)
-        b = self.np.conj(a)
+        a = (1 + 1j) / 2
+        b = (1 - 1j) / 2
         return self._cast(
             [
-                [1, 0, 0, 0],
-                [0, 1, 0, 0],
-                [0, 0, a, b],
+                [1 + 0j, 0, 0, 0],
+                [0, 1 + 0j, 0, 0],
                 [0, 0, b, a],
+                [0, 0, a, b],
             ],
             dtype=self.dtype,
         )
 
     def CRX(self, theta):
-        theta = self._cast_parameter(theta)
         cos = self.np.cos(theta / 2.0) + 0j
         isin = -1j * self.np.sin(theta / 2.0)
         matrix = [
@@ -218,14 +198,12 @@ def CRX(self, theta):
         return self._cast(matrix, dtype=self.dtype)
 
     def CRY(self, theta):
-        theta = self._cast_parameter(theta)
         cos = self.np.cos(theta / 2.0) + 0j
         sin = self.np.sin(theta / 2.0) + 0j
         matrix = [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, cos, -sin], [0, 0, sin, cos]]
         return self._cast(matrix, dtype=self.dtype)
 
     def CRZ(self, theta):
-        theta = self._cast_parameter(theta)
         phase = self.np.exp(0.5j * theta)
         matrix = [
             [1, 0, 0, 0],
@@ -236,7 +214,6 @@ def CRZ(self, theta):
         return self._cast(matrix, dtype=self.dtype)
 
     def CU1(self, theta):
-        theta = self._cast_parameter(theta)
         phase = self.np.exp(1j * theta)
         matrix = [
             [1, 0, 0, 0],
@@ -247,8 +224,6 @@ def CU1(self, theta):
         return self._cast(matrix, dtype=self.dtype)
 
     def CU2(self, phi, lam):
-        phi = self._cast_parameter(phi)
-        lam = self._cast_parameter(lam)
         eplus = self.np.exp(1j * (phi + lam) / 2.0) / math.sqrt(2)
         eminus = self.np.exp(1j * (phi - lam) / 2.0) / math.sqrt(2)
         matrix = [
@@ -260,9 +235,6 @@ def CU2(self, phi, lam):
         return self._cast(matrix, dtype=self.dtype)
 
     def CU3(self, theta, phi, lam):
-        theta = self._cast_parameter(theta)
-        phi = self._cast_parameter(phi)
-        lam = self._cast_parameter(lam)
         cost = self.np.cos(theta / 2)
         sint = self.np.sin(theta / 2)
         eplus = self.np.exp(1j * (phi + lam) / 2.0)
@@ -324,8 +296,6 @@ def FSWAP(self):
         )
 
     def fSim(self, theta, phi):
-        theta = self._cast_parameter(theta)
-        phi = self._cast_parameter(phi)
         cost = self.np.cos(theta) + 0j
         isint = -1j * self.np.sin(theta)
         phase = self.np.exp(-1j * phi)
@@ -355,7 +325,6 @@ def SYC(self):
         )
 
     def GeneralizedfSim(self, u, phi):
-        phi = self._cast_parameter(phi)
         phase = self.np.exp(-1j * phi)
         return self._cast(
             [
@@ -368,7 +337,6 @@ def GeneralizedfSim(self, u, phi):
         )
 
     def RXX(self, theta):
-        theta = self._cast_parameter(theta)
         cos = self.np.cos(theta / 2.0) + 0j
         isin = -1j * self.np.sin(theta / 2.0)
         return self._cast(
@@ -382,7 +350,6 @@ def RXX(self, theta):
         )
 
     def RYY(self, theta):
-        theta = self._cast_parameter(theta)
         cos = self.np.cos(theta / 2.0) + 0j
         isin = -1j * self.np.sin(theta / 2.0)
         return self._cast(
@@ -396,7 +363,6 @@ def RYY(self, theta):
         )
 
     def RZZ(self, theta):
-        theta = self._cast_parameter(theta)
         phase = self.np.exp(0.5j * theta)
         return self._cast(
             [
@@ -409,7 +375,6 @@ def RZZ(self, theta):
         )
 
     def RZX(self, theta):
-        theta = self._cast_parameter(theta)
         cos, sin = self.np.cos(theta / 2) + 0j, self.np.sin(theta / 2) + 0j
         return self._cast(
             [
@@ -422,7 +387,6 @@ def RZX(self, theta):
         )
 
     def RXXYY(self, theta):
-        theta = self._cast_parameter(theta)
         cos, sin = self.np.cos(theta / 2) + 0j, self.np.sin(theta / 2) + 0j
         return self._cast(
             [
@@ -435,11 +399,6 @@ def RXXYY(self, theta):
         )
 
     def MS(self, phi0, phi1, theta):
-        phi0, phi1, theta = (
-            self._cast_parameter(phi0),
-            self._cast_parameter(phi1),
-            self._cast_parameter(theta),
-        )
         plus = self.np.exp(1.0j * (phi0 + phi1))
         minus = self.np.exp(1.0j * (phi0 - phi1))
         cos = self.np.cos(theta / 2) + 0j
@@ -455,7 +414,6 @@ def MS(self, phi0, phi1, theta):
         )
 
     def GIVENS(self, theta):
-        theta = self._cast_parameter(theta)
         return self._cast(
             [
                 [1, 0, 0, 0],
@@ -514,7 +472,6 @@ def CCZ(self):
         )
 
     def DEUTSCH(self, theta):
-        theta = self._cast_parameter(theta)
         sin = self.np.sin(theta) + 0j  # 0j necessary for right tensorflow dtype
         cos = self.np.cos(theta) + 0j
         return self._cast(
@@ -532,8 +489,6 @@ def DEUTSCH(self, theta):
         )
 
     def GeneralizedRBS(self, qubits_in, qubits_out, theta, phi):
-        theta = self._cast_parameter(theta)
-        phi = self._cast_parameter(phi)
         bitstring_length = len(qubits_in) + len(qubits_out)
         integer_in = "".join(
             ["1" if k in qubits_in else "0" for k in range(bitstring_length)]
diff --git a/src/qibo/backends/numpy.py b/src/qibo/backends/numpy.py
index 21be02bfce..bda83e75cd 100644
--- a/src/qibo/backends/numpy.py
+++ b/src/qibo/backends/numpy.py
@@ -433,7 +433,6 @@ def execute_circuit(self, circuit, initial_state=None, nshots=1000):
                 if initial_state is None:
                     state = self.zero_density_matrix(nqubits)
                 else:
-                    # cast to proper complex type
                     state = self.cast(initial_state)
 
                 for gate in circuit.queue:
@@ -443,7 +442,6 @@ def execute_circuit(self, circuit, initial_state=None, nshots=1000):
                 if initial_state is None:
                     state = self.zero_state(nqubits)
                 else:
-                    # cast to proper complex type
                     state = self.cast(initial_state)
 
                 for gate in circuit.queue:
diff --git a/src/qibo/backends/pytorch.py b/src/qibo/backends/pytorch.py
index 186c61d600..ed4da5bfbd 100644
--- a/src/qibo/backends/pytorch.py
+++ b/src/qibo/backends/pytorch.py
@@ -14,25 +14,20 @@ class TorchMatrices(NumpyMatrices):
 
     Args:
         dtype (torch.dtype): Data type of the matrices.
-        requires_grad (bool): If ``True`` the matrices require gradient.
     """
 
-    def __init__(self, dtype, requires_grad):
+    def __init__(self, dtype):
         import torch  # pylint: disable=import-outside-toplevel  # type: ignore
 
         super().__init__(dtype)
         self.np = torch
         self.dtype = dtype
-        self.requires_grad = requires_grad
 
     def _cast(self, x, dtype):
         flattened = [item for sublist in x for item in sublist]
         tensor_list = [self.np.as_tensor(i, dtype=dtype) for i in flattened]
         return self.np.stack(tensor_list).reshape(len(x), len(x))
 
-    def _cast_parameter(self, x):
-        return self.np.tensor(x, dtype=self.dtype, requires_grad=self.requires_grad)
-
     def Unitary(self, u):
         return self._cast(u, dtype=self.dtype)
 
@@ -42,9 +37,6 @@ def __init__(self):
         super().__init__()
         import torch  # pylint: disable=import-outside-toplevel  # type: ignore
 
-        # Global variable to enable or disable gradient calculation
-        self.gradients = True
-
         self.np = torch
 
         self.name = "pytorch"
@@ -54,8 +46,11 @@ def __init__(self):
             "torch": self.np.__version__,
         }
 
+        # Default data type used for the gate matrices is complex128
         self.dtype = self._torch_dtype(self.dtype)
-        self.matrices = TorchMatrices(self.dtype, requires_grad=self.gradients)
+        # Default data type used for the real gate parameters is float64
+        self.parameter_dtype = self._torch_dtype("float64")
+        self.matrices = TorchMatrices(self.dtype)
         self.device = self.np.device("cuda:0" if torch.cuda.is_available() else "cpu")
         self.nthreads = 0
         self.tensor_types = (self.np.Tensor, np.ndarray)
@@ -63,17 +58,13 @@ def __init__(self):
         # These functions in Torch works in a different way than numpy or have different names
         self.np.transpose = self.np.permute
         self.np.copy = self.np.clone
+        self.np.power = self.np.pow
         self.np.expand_dims = self.np.unsqueeze
         self.np.mod = self.np.remainder
         self.np.right_shift = self.np.bitwise_right_shift
         self.np.sign = self.np.sgn
         self.np.flatnonzero = lambda x: self.np.nonzero(x).flatten()
 
-    def requires_grad(self, requires_grad):
-        """Enable or disable gradient calculation."""
-        self.gradients = requires_grad
-        self.matrices.requires_grad = requires_grad
-
     def _torch_dtype(self, dtype):
         if dtype == "float":
             dtype += "32"
@@ -87,7 +78,6 @@ def cast(
         x,
         dtype=None,
         copy: bool = False,
-        requires_grad: bool = None,
     ):
         """Casts input as a Torch tensor of the specified dtype.
 
@@ -102,12 +92,7 @@ def cast(
                 Defaults to ``None``.
             copy (bool, optional): If ``True``, the input tensor is copied before casting.
                 Defaults to ``False``.
-            requires_grad (bool, optional): If ``True``, the input tensor requires gradient.
-                If ``False``, the input tensor does not require gradient.
-                If ``None``, the default gradient setting of the backend is used.
         """
-        if requires_grad is None:
-            requires_grad = self.gradients
 
         if dtype is None:
             dtype = self.dtype
@@ -116,10 +101,6 @@ def cast(
         elif not isinstance(dtype, self.np.dtype):
             dtype = self._torch_dtype(str(dtype))
 
-        # check if dtype is an integer to remove gradients
-        if dtype in [self.np.int32, self.np.int64, self.np.int8, self.np.int16]:
-            requires_grad = False
-
         if isinstance(x, self.np.Tensor):
             x = x.to(dtype)
         elif (
@@ -129,13 +110,59 @@ def cast(
         ):
             x = self.np.stack(x)
         else:
-            x = self.np.tensor(x, dtype=dtype, requires_grad=requires_grad)
+            x = self.np.tensor(x, dtype=dtype)
 
         if copy:
             return x.clone()
-
         return x
 
+    def matrix_parametrized(self, gate):
+        """Convert a parametrized gate to its matrix representation in the computational basis."""
+        name = gate.__class__.__name__
+        _matrix = getattr(self.matrices, name)
+        if name == "GeneralizedRBS":
+            for parameter in ["theta", "phi"]:
+                if not isinstance(gate.init_kwargs[parameter], self.np.Tensor):
+                    gate.init_kwargs[parameter] = self._cast_parameter(
+                        gate.init_kwargs[parameter], trainable=gate.trainable
+                    )
+
+            _matrix = _matrix(
+                qubits_in=gate.init_args[0],
+                qubits_out=gate.init_args[1],
+                theta=gate.init_kwargs["theta"],
+                phi=gate.init_kwargs["phi"],
+            )
+            return _matrix
+        else:
+            new_parameters = []
+            for parameter in gate.parameters:
+                if not isinstance(parameter, self.np.Tensor):
+                    parameter = self._cast_parameter(
+                        parameter, trainable=gate.trainable
+                    )
+                elif parameter.requires_grad:
+                    gate.trainable = True
+                new_parameters.append(parameter)
+            gate.parameters = tuple(new_parameters)
+        _matrix = _matrix(*gate.parameters)
+        return _matrix
+
+    def _cast_parameter(self, x, trainable):
+        """Cast a gate parameter to a torch tensor.
+
+        Args:
+            x (Union[int, float, complex]): Parameter to be casted.
+            trainable (bool): If ``True``, the tensor requires gradient.
+        """
+        if isinstance(x, int) and trainable:
+            return self.np.tensor(x, dtype=self.parameter_dtype, requires_grad=True)
+        if isinstance(x, float):
+            return self.np.tensor(
+                x, dtype=self.parameter_dtype, requires_grad=trainable
+            )
+        return self.np.tensor(x, dtype=self.dtype, requires_grad=trainable)
+
     def is_sparse(self, x):
         if isinstance(x, self.np.Tensor):
             return x.is_sparse
diff --git a/src/qibo/gates/channels.py b/src/qibo/gates/channels.py
index be43678318..ba4f285df2 100644
--- a/src/qibo/gates/channels.py
+++ b/src/qibo/gates/channels.py
@@ -653,12 +653,14 @@ def apply_density_matrix(self, backend, state, nqubits):
                 self.init_kwargs["p_1"],
                 self.init_kwargs["e_t2"],
             )
-            matrix = [
-                [1 - preset1, 0, 0, preset0],
-                [0, e_t2, 0, 0],
-                [0, 0, e_t2, 0],
-                [preset1, 0, 0, 1 - preset0],
-            ]
+            matrix = backend.cast(
+                [
+                    [1 - preset1, 0, 0, preset0],
+                    [0, e_t2, 0, 0],
+                    [0, 0, e_t2, 0],
+                    [preset1, 0, 0, 1 - preset0],
+                ]
+            )
 
             qubits = (qubit, qubit + nqubits)
             gate = Unitary(matrix, *qubits)
diff --git a/src/qibo/gates/gates.py b/src/qibo/gates/gates.py
index 98d6baa289..618279c897 100644
--- a/src/qibo/gates/gates.py
+++ b/src/qibo/gates/gates.py
@@ -506,13 +506,16 @@ def qasm_label(self):
 class Align(ParametrizedGate):
     """Aligns proceeding qubit operations and (optionally) waits ``delay`` amount of time.
 
+    .. note::
+        For this gate, the ``trainable`` parameter is by default set to ``False``.
+
     Args:
         q (int): The qubit ID.
         delay (int, optional): The time (in ns) for which to delay circuit execution on the specified qubits.
             Defaults to ``0`` (zero).
     """
 
-    def __init__(self, q, delay=0, trainable=True):
+    def __init__(self, q, delay=0, trainable=False):
         if not isinstance(delay, int):
             raise_error(
                 TypeError, f"delay must be type int, but it is type {type(delay)}."
@@ -2600,9 +2603,10 @@ def __init__(
     @Gate.parameters.setter
     def parameters(self, x):
         shape = self.parameters[0].shape
-        engine = _check_engine(self.parameters[0])
+        engine = _check_engine(x)
         # Reshape doesn't accept a tuple if engine is pytorch.
-        x = x[0] if type(x) is tuple else x
+        if isinstance(x, tuple):
+            x = x[0]
         self._parameters = (engine.reshape(x, shape),)
         for gate in self.device_gates:  # pragma: no cover
             gate.parameters = x
@@ -2633,7 +2637,9 @@ def _dagger(self):
 
 def _check_engine(array):
     """Check if the array is a numpy or torch tensor and return the corresponding library."""
-    if array.__class__.__name__ == "Tensor":
+    if (array.__class__.__name__ == "Tensor") or (
+        isinstance(array, tuple) and array[0].__class__.__name__ == "Tensor"
+    ):
         import torch  # pylint: disable=C0415
 
         return torch
diff --git a/src/qibo/models/circuit.py b/src/qibo/models/circuit.py
index 8bfb423abf..35d18100a6 100644
--- a/src/qibo/models/circuit.py
+++ b/src/qibo/models/circuit.py
@@ -1142,6 +1142,9 @@ def from_dict(cls, raw):
     def to_qasm(self):
         """Convert circuit to QASM.
 
+        .. note::
+            This method does not support multi-controlled gates and gates with ``torch.Tensor`` as parameters.
+
         Args:
             filename (str): The filename where the code is saved.
         """
@@ -1174,7 +1177,7 @@ def to_qasm(self):
 
             qubits = ",".join(f"q[{i}]" for i in gate.qubits)
             if isinstance(gate, gates.ParametrizedGate):
-                params = (str(x) for x in gate.parameters)
+                params = (str(float(x)) for x in gate.parameters)
                 name = f"{gate.qasm_label}({', '.join(params)})"
             else:
                 name = gate.qasm_label
diff --git a/src/qibo/models/error_mitigation.py b/src/qibo/models/error_mitigation.py
index e85b7cfd73..60a7da94af 100644
--- a/src/qibo/models/error_mitigation.py
+++ b/src/qibo/models/error_mitigation.py
@@ -330,6 +330,7 @@ def _curve_fit(
     if backend.name == "pytorch":
         # pytorch has some problems with the `scipy.optim.curve_fit` function
         # thus we use a `torch.optim` optimizer
+        params.requires_grad = True
         loss = lambda pred, target: backend.np.mean((pred - target) ** 2)
         optimizer = backend.np.optim.LBFGS(
             [params], lr=lr, max_iter=max_iter, tolerance_grad=tolerance_grad
diff --git a/src/qibo/transpiler/decompositions.py b/src/qibo/transpiler/decompositions.py
index 29d691402a..693c815e29 100644
--- a/src/qibo/transpiler/decompositions.py
+++ b/src/qibo/transpiler/decompositions.py
@@ -1,14 +1,11 @@
 import numpy as np
 
 from qibo import gates
-from qibo.backends import NumpyBackend
 from qibo.transpiler.unitary_decompositions import (
     two_qubit_decomposition,
     u3_decomposition,
 )
 
-backend = NumpyBackend()
-
 
 class GateDecompositions:
     """Abstract data structure that holds decompositions of gates."""
@@ -20,27 +17,35 @@ def add(self, gate, decomposition):
         """Register a decomposition for a gate."""
         self.decompositions[gate] = decomposition
 
-    def count_2q(self, gate):
-        """Count the number of two-qubit gates in the decomposition of the given gate."""
+    def _check_instance(self, gate, backend=None):
+        special_gates = (
+            gates.FusedGate,
+            gates.Unitary,
+            gates.GeneralizedfSim,
+            gates.fSim,
+        )
+        decomposition = self.decompositions[gate.__class__]
         if gate.parameters:
-            decomposition = self.decompositions[gate.__class__](gate)
-        else:
-            decomposition = self.decompositions[gate.__class__]
+            decomposition = (
+                decomposition(gate, backend)
+                if isinstance(gate, special_gates)
+                else decomposition(gate)
+            )
+        return decomposition
+
+    def count_2q(self, gate, backend):
+        """Count the number of two-qubit gates in the decomposition of the given gate."""
+        decomposition = self._check_instance(gate, backend)
         return len(tuple(g for g in decomposition if len(g.qubits) > 1))
 
-    def count_1q(self, gate):
+    def count_1q(self, gate, backend):
         """Count the number of single qubit gates in the decomposition of the given gate."""
-        if gate.parameters:
-            decomposition = self.decompositions[gate.__class__](gate)
-        else:
-            decomposition = self.decompositions[gate.__class__]
+        decomposition = self._check_instance(gate, backend)
         return len(tuple(g for g in decomposition if len(g.qubits) == 1))
 
-    def __call__(self, gate):
+    def __call__(self, gate, backend=None):
         """Decompose a gate."""
-        decomposition = self.decompositions[gate.__class__]
-        if callable(decomposition):
-            decomposition = decomposition(gate)
+        decomposition = self._check_instance(gate, backend)
         return [
             g.on_qubits({i: q for i, q in enumerate(gate.qubits)})
             for g in decomposition
@@ -48,15 +53,17 @@ def __call__(self, gate):
 
 
 def _u3_to_gpi2(t, p, l):
-    """Decompose a U3 gate into GPI2 gates, the decomposition is optimized to use the minimum number of gates..
+    """Decompose a :class:`qibo.gates.U3` gate into :class:`qibo.gates.GPI2` gates.
+
+    The decomposition is optimized to use the minimum number of gates.
 
     Args:
-        t (float): theta parameter of U3 gate.
-        p (float): phi parameter of U3 gate.
-        l (float): lambda parameter of U3 gate.
+        t (float): first parameter of :class:`qibo.gates.U3` gate.
+        p (float): second parameter of :class:`qibo.gates.U3` gate.
+        l (float): third parameter of :class:`qibo.gates.U3` gate.
 
     Returns:
-        decomposition (list): list of native gates that decompose the U3 gate.
+        list: Native gates that decompose the :class:`qibo.gates.U3` gate.
     """
     decomposition = []
     if l != 0.0:
@@ -108,10 +115,12 @@ def _u3_to_gpi2(t, p, l):
 )
 gpi2_dec.add(gates.U3, lambda gate: _u3_to_gpi2(*gate.parameters))
 gpi2_dec.add(
-    gates.Unitary, lambda gate: _u3_to_gpi2(*u3_decomposition(gate.parameters[0]))
+    gates.Unitary,
+    lambda gate, backend: _u3_to_gpi2(*u3_decomposition(gate.parameters[0], backend)),
 )
 gpi2_dec.add(
-    gates.FusedGate, lambda gate: _u3_to_gpi2(*u3_decomposition(gate.matrix(backend)))
+    gates.FusedGate,
+    lambda gate, backend: _u3_to_gpi2(*u3_decomposition(gate.matrix(backend), backend)),
 )
 
 # Decompose single qubit gates using U3
@@ -126,7 +135,8 @@ def _u3_to_gpi2(t, p, l):
 u3_dec.add(gates.TDG, [gates.RZ(0, -np.pi / 4)])
 u3_dec.add(gates.SX, [gates.U3(0, np.pi / 2, -np.pi / 2, np.pi / 2)])
 u3_dec.add(
-    gates.RX, lambda gate: [gates.U3(0, gate.parameters[0], -np.pi / 2, np.pi / 2)]
+    gates.RX,
+    lambda gate: [gates.U3(0, gate.parameters[0], -np.pi / 2, np.pi / 2)],
 )
 u3_dec.add(gates.RY, lambda gate: [gates.U3(0, gate.parameters[0], 0, 0)])
 u3_dec.add(gates.RZ, lambda gate: [gates.RZ(0, gate.parameters[0])])
@@ -139,7 +149,12 @@ def _u3_to_gpi2(t, p, l):
     ],
 )
 u3_dec.add(
-    gates.GPI2, lambda gate: [gates.U3(0, *u3_decomposition(gate.matrix(backend)))]
+    gates.GPI2,
+    lambda gate: [
+        gates.U3(
+            0, np.pi / 2, gate.parameters[0] - np.pi / 2, np.pi / 2 - gate.parameters[0]
+        ),
+    ],
 )
 u3_dec.add(gates.U1, lambda gate: [gates.RZ(0, gate.parameters[0])])
 u3_dec.add(
@@ -154,11 +169,13 @@ def _u3_to_gpi2(t, p, l):
 )
 u3_dec.add(
     gates.Unitary,
-    lambda gate: [gates.U3(0, *u3_decomposition(gate.parameters[0]))],
+    lambda gate, backend: [gates.U3(0, *u3_decomposition(gate.parameters[0], backend))],
 )
 u3_dec.add(
     gates.FusedGate,
-    lambda gate: [gates.U3(0, *u3_decomposition(gate.matrix(backend)))],
+    lambda gate, backend: [
+        gates.U3(0, *u3_decomposition(gate.matrix(backend), backend))
+    ],
 )
 
 # register the iSWAP decompositions
@@ -376,15 +393,21 @@ def _u3_to_gpi2(t, p, l):
 )
 cz_dec.add(
     gates.Unitary,
-    lambda gate: two_qubit_decomposition(0, 1, gate.parameters[0], backend=backend),
+    lambda gate, backend: two_qubit_decomposition(
+        0, 1, gate.parameters[0], backend=backend
+    ),
 )
 cz_dec.add(
     gates.fSim,
-    lambda gate: two_qubit_decomposition(0, 1, gate.matrix(backend), backend=backend),
+    lambda gate, backend: two_qubit_decomposition(
+        0, 1, gate.matrix(backend), backend=backend
+    ),
 )
 cz_dec.add(
     gates.GeneralizedfSim,
-    lambda gate: two_qubit_decomposition(0, 1, gate.matrix(backend), backend=backend),
+    lambda gate, backend: two_qubit_decomposition(
+        0, 1, gate.matrix(backend), backend=backend
+    ),
 )
 
 # temporary CNOT decompositions for CNOT, CZ, SWAP
diff --git a/src/qibo/transpiler/unitary_decompositions.py b/src/qibo/transpiler/unitary_decompositions.py
index ad892db734..74fb7f2e02 100644
--- a/src/qibo/transpiler/unitary_decompositions.py
+++ b/src/qibo/transpiler/unitary_decompositions.py
@@ -1,8 +1,7 @@
 import numpy as np
 
 from qibo import gates, matrices
-from qibo.backends import _check_backend
-from qibo.config import PRECISION_TOL, raise_error
+from qibo.config import raise_error
 
 magic_basis = np.array(
     [[1, -1j, 0, 0], [0, 0, 1, -1j], [0, 0, -1, -1j], [1, 1j, 0, 0]]
@@ -15,7 +14,7 @@
 H = np.array([[1, 1], [1, -1]]) / np.sqrt(2)
 
 
-def u3_decomposition(unitary):
+def u3_decomposition(unitary, backend):
     """Decomposes arbitrary one-qubit gates to U3.
 
     Args:
@@ -24,18 +23,19 @@ def u3_decomposition(unitary):
     Returns:
         (float, float, float): parameters of U3 gate.
     """
+    unitary = backend.cast(unitary)
     # https://github.com/Qiskit/qiskit-terra/blob/d2e3340adb79719f9154b665e8f6d8dc26b3e0aa/qiskit/quantum_info/synthesis/one_qubit_decompose.py#L221
-    su2 = unitary / np.sqrt(np.linalg.det(unitary))
-    theta = 2 * np.arctan2(abs(su2[1, 0]), abs(su2[0, 0]))
-    plus = np.angle(su2[1, 1])
-    minus = np.angle(su2[1, 0])
+    su2 = unitary / backend.np.sqrt(backend.np.linalg.det(unitary))
+    theta = 2 * backend.np.arctan2(backend.np.abs(su2[1, 0]), backend.np.abs(su2[0, 0]))
+    plus = backend.np.angle(su2[1, 1])
+    minus = backend.np.angle(su2[1, 0])
     phi = plus + minus
     lam = plus - minus
 
     return theta, phi, lam
 
 
-def calculate_psi(unitary, magic_basis=magic_basis, backend=None):
+def calculate_psi(unitary, backend, magic_basis=magic_basis):
     """Solves the eigenvalue problem of :math:`U^{T} U`.
 
     See step (1) of Appendix A in arXiv:quant-ph/0011050.
@@ -50,7 +50,6 @@ def calculate_psi(unitary, magic_basis=magic_basis, backend=None):
     Returns:
         ndarray: Eigenvectors in the computational basis and eigenvalues of :math:`U^{T} U`.
     """
-    backend = _check_backend(backend)
 
     if backend.__class__.__name__ in [
         "CupyBackend",
@@ -84,7 +83,7 @@ def calculate_psi(unitary, magic_basis=magic_basis, backend=None):
     return psi, eigvals
 
 
-def schmidt_decompose(state, backend=None):
+def schmidt_decompose(state, backend):
     """Decomposes a two-qubit product state to its single-qubit parts.
 
     Args:
@@ -94,7 +93,6 @@ def schmidt_decompose(state, backend=None):
         (ndarray, ndarray): decomposition
 
     """
-    backend = _check_backend(backend)
     # tf.linalg.svd has a different behaviour
     if backend.__class__.__name__ == "TensorflowBackend":
         u, d, v = np.linalg.svd(backend.np.reshape(state, (2, 2)))
@@ -108,7 +106,7 @@ def schmidt_decompose(state, backend=None):
     return u[:, 0], v[0]
 
 
-def calculate_single_qubit_unitaries(psi, backend=None):
+def calculate_single_qubit_unitaries(psi, backend):
     """Calculates local unitaries that maps a maximally entangled basis to the magic basis.
 
     See Lemma 1 of Appendix A in arXiv:quant-ph/0011050.
@@ -119,13 +117,12 @@ def calculate_single_qubit_unitaries(psi, backend=None):
     Returns:
         (ndarray, ndarray): Local unitaries UA and UB that map the given basis to the magic basis.
     """
-    backend = _check_backend(backend)
     psi_magic = backend.np.matmul(backend.np.conj(backend.cast(magic_basis)).T, psi)
     if (
         backend.np.real(
             backend.calculate_norm_density_matrix(backend.np.imag(psi_magic))
         )
-        > PRECISION_TOL
+        > 1e-6
     ):  # pragma: no cover
         raise_error(NotImplementedError, "Given state is not real in the magic basis.")
     psi_bar = backend.cast(psi.T, copy=True)
@@ -154,13 +151,12 @@ def calculate_single_qubit_unitaries(psi, backend=None):
     return ua, ub
 
 
-def calculate_diagonal(unitary, ua, ub, va, vb, backend=None):
+def calculate_diagonal(unitary, ua, ub, va, vb, backend):
     """Calculates Ud matrix that can be written as exp(-iH).
 
     See Eq. (A1) in arXiv:quant-ph/0011050.
     Ud is diagonal in the magic and Bell basis.
     """
-    backend = _check_backend(backend)
     # normalize U_A, U_B, V_A, V_B so that detU_d = 1
     # this is required so that sum(lambdas) = 0
     # and Ud can be written as exp(-iH)
@@ -188,7 +184,7 @@ def calculate_diagonal(unitary, ua, ub, va, vb, backend=None):
 
 def magic_decomposition(unitary, backend=None):
     """Decomposes an arbitrary unitary to (A1) from arXiv:quant-ph/0011050."""
-    backend = _check_backend(backend)
+
     unitary = backend.cast(unitary)
     psi, eigvals = calculate_psi(unitary, backend=backend)
     psi_tilde = backend.np.conj(backend.np.sqrt(eigvals)) * backend.np.matmul(
@@ -202,10 +198,8 @@ def magic_decomposition(unitary, backend=None):
     return calculate_diagonal(unitary, ua, ub, va, vb, backend=backend)
 
 
-def to_bell_diagonal(ud, bell_basis=bell_basis, backend=None):
+def to_bell_diagonal(ud, backend, bell_basis=bell_basis):
     """Transforms a matrix to the Bell basis and checks if it is diagonal."""
-    backend = _check_backend(backend)
-
     ud = backend.cast(ud)
     bell_basis = backend.cast(bell_basis)
 
@@ -213,20 +207,21 @@ def to_bell_diagonal(ud, bell_basis=bell_basis, backend=None):
         backend.np.transpose(backend.np.conj(bell_basis), (1, 0)) @ ud @ bell_basis
     )
     ud_diag = backend.np.diag(ud_bell)
-    if not backend.np.allclose(backend.np.diag(ud_diag), ud_bell):  # pragma: no cover
+    if not backend.np.allclose(
+        backend.np.diag(ud_diag), ud_bell, atol=1e-6, rtol=1e-6
+    ):  # pragma: no cover
         return None
-    uprod = backend.np.prod(ud_diag)
-    if not np.allclose(backend.to_numpy(uprod), 1):  # pragma: no cover
+    uprod = backend.to_numpy(backend.np.prod(ud_diag))
+    if not np.allclose(uprod, 1.0, atol=1e-6, rtol=1e-6):  # pragma: no cover
         return None
     return ud_diag
 
 
-def calculate_h_vector(ud_diag, backend=None):
+def calculate_h_vector(ud_diag, backend):
     """Finds h parameters corresponding to exp(-iH).
 
     See Eq. (4)-(5) in arXiv:quant-ph/0307177.
     """
-    backend = _check_backend(backend)
     lambdas = -backend.np.angle(ud_diag)
     hx = (lambdas[0] + lambdas[2]) / 2.0
     hy = (lambdas[1] + lambdas[2]) / 2.0
@@ -234,9 +229,8 @@ def calculate_h_vector(ud_diag, backend=None):
     return hx, hy, hz
 
 
-def cnot_decomposition(q0, q1, hx, hy, hz, backend=None):
+def cnot_decomposition(q0, q1, hx, hy, hz, backend):
     """Performs decomposition (6) from arXiv:quant-ph/0307177."""
-    backend = _check_backend(backend)
     h = backend.cast(H)
     u3 = backend.cast(-1j * matrices.H)
     # use corrected version from PRA paper (not arXiv)
@@ -260,9 +254,8 @@ def cnot_decomposition(q0, q1, hx, hy, hz, backend=None):
     ]
 
 
-def cnot_decomposition_light(q0, q1, hx, hy, backend=None):
+def cnot_decomposition_light(q0, q1, hx, hy, backend):
     """Performs decomposition (24) from arXiv:quant-ph/0307177."""
-    backend = _check_backend(backend)
     h = backend.cast(H)
     w = backend.cast((matrices.I - 1j * matrices.X) / np.sqrt(2))
     u2 = gates.RX(0, 2 * hx).matrix(backend)
@@ -280,18 +273,18 @@ def cnot_decomposition_light(q0, q1, hx, hy, backend=None):
     ]
 
 
-def two_qubit_decomposition(q0, q1, unitary, backend=None):
+def two_qubit_decomposition(q0, q1, unitary, backend):
     """Performs two qubit unitary gate decomposition (24) from arXiv:quant-ph/0307177.
 
     Args:
         q0 (int): index of the first qubit.
         q1 (int): index of the second qubit.
         unitary (ndarray): Unitary :math:`4 \\times 4` to be decomposed.
+        backend (:class:`qibo.backends.Backend`): Backend to use for calculations.
 
     Returns:
         (list): gates implementing decomposition (24) from arXiv:quant-ph/0307177
     """
-    backend = _check_backend(backend)
 
     ud_diag = to_bell_diagonal(unitary, backend=backend)
     ud = None
diff --git a/src/qibo/transpiler/unroller.py b/src/qibo/transpiler/unroller.py
index 9032e5e306..d3b564884b 100644
--- a/src/qibo/transpiler/unroller.py
+++ b/src/qibo/transpiler/unroller.py
@@ -1,6 +1,7 @@
 from enum import Flag, auto
 
 from qibo import gates
+from qibo.backends import _check_backend
 from qibo.config import raise_error
 from qibo.models import Circuit
 from qibo.transpiler._exceptions import DecompositionError
@@ -75,6 +76,7 @@ class Unroller:
     Args:
         native_gates (:class:`qibo.transpiler.unroller.NativeGates`): native gates to use
             in the transpiled circuit.
+        backend (:class:`qibo.backends.Backend`): backend to use for gate matrix.
 
     Returns:
         (:class:`qibo.models.circuit.Circuit`): equivalent circuit with native gates.
@@ -83,8 +85,10 @@ class Unroller:
     def __init__(
         self,
         native_gates: NativeGates,
+        backend=None,
     ):
         self.native_gates = native_gates
+        self.backend = backend
 
     def __call__(self, circuit: Circuit):
         """Decomposes a circuit into native gates.
@@ -101,6 +105,7 @@ def __call__(self, circuit: Circuit):
                 translate_gate(
                     gate,
                     self.native_gates,
+                    backend=self.backend,
                 )
             )
         return translated_circuit
@@ -142,6 +147,7 @@ def assert_decomposition(
 def translate_gate(
     gate,
     native_gates: NativeGates,
+    backend=None,
 ):
     """Maps gates to a hardware-native implementation.
 
@@ -149,10 +155,13 @@ def translate_gate(
         gate (:class:`qibo.gates.abstract.Gate`): gate to be decomposed.
         native_gates (:class:`qibo.transpiler.unroller.NativeGates`): native gates
             to use in the decomposition.
+        backend (:class:`qibo.backends.Backend`): backend to use for gate matrix.
 
     Returns:
         list: List of native gates that decompose the input gate.
     """
+    backend = _check_backend(backend)
+
     if isinstance(gate, (gates.I, gates.Align)):
         return gate
 
@@ -162,21 +171,23 @@ def translate_gate(
         return gate
 
     if len(gate.qubits) == 1:
-        return _translate_single_qubit_gates(gate, native_gates)
+        return _translate_single_qubit_gates(gate, native_gates, backend)
 
-    decomposition_2q = _translate_two_qubit_gates(gate, native_gates)
+    decomposition_2q = _translate_two_qubit_gates(gate, native_gates, backend)
     final_decomposition = []
     for decomposed_2q_gate in decomposition_2q:
         if len(decomposed_2q_gate.qubits) == 1:
             final_decomposition += _translate_single_qubit_gates(
-                decomposed_2q_gate, native_gates
+                decomposed_2q_gate, native_gates, backend
             )
         else:
             final_decomposition.append(decomposed_2q_gate)
     return final_decomposition
 
 
-def _translate_single_qubit_gates(gate: gates.Gate, single_qubit_natives: NativeGates):
+def _translate_single_qubit_gates(
+    gate: gates.Gate, single_qubit_natives: NativeGates, backend
+):
     """Helper method for :meth:`translate_gate`.
 
     Maps single qubit gates to a hardware-native implementation.
@@ -185,20 +196,21 @@ def _translate_single_qubit_gates(gate: gates.Gate, single_qubit_natives: Native
         gate (:class:`qibo.gates.abstract.Gate`): gate to be decomposed.
         single_qubit_natives (:class:`qibo.transpiler.unroller.NativeGates`): single
             qubit native gates.
+        backend (:class:`qibo.backends.Backend`): backend to use for gate matrix.
 
     Returns:
         list: List of native gates that decompose the input gate.
     """
     if NativeGates.U3 & single_qubit_natives:
-        return u3_dec(gate)
+        return u3_dec(gate, backend)
 
     if NativeGates.GPI2 & single_qubit_natives:
-        return gpi2_dec(gate)
+        return gpi2_dec(gate, backend)
 
     raise_error(DecompositionError, "Use U3 or GPI2 as single qubit native gates")
 
 
-def _translate_two_qubit_gates(gate: gates.Gate, native_gates: NativeGates):
+def _translate_two_qubit_gates(gate: gates.Gate, native_gates: NativeGates, backend):
     """Helper method for :meth:`translate_gate`.
 
     Maps two qubit gates to a hardware-native implementation.
@@ -207,6 +219,7 @@ def _translate_two_qubit_gates(gate: gates.Gate, native_gates: NativeGates):
         gate (:class:`qibo.gates.abstract.Gate`): gate to be decomposed.
         native_gates (:class:`qibo.transpiler.unroller.NativeGates`): native gates
             supported by the quantum hardware.
+        backend (:class:`qibo.backends.Backend`): backend to use for gate matrix.
 
     Returns:
         list: List of native gates that decompose the input gate.
@@ -216,36 +229,40 @@ def _translate_two_qubit_gates(gate: gates.Gate, native_gates: NativeGates):
     ) is NativeGates.CZ | NativeGates.iSWAP:
         # Check for a special optimized decomposition.
         if gate.__class__ in opt_dec.decompositions:
-            return opt_dec(gate)
+            return opt_dec(gate, backend)
         # Check if the gate has a CZ decomposition
         if not gate.__class__ in iswap_dec.decompositions:
-            return cz_dec(gate)
+            return cz_dec(gate, backend)
         # Check the decomposition with less 2 qubit gates.
 
-        if cz_dec.count_2q(gate) < iswap_dec.count_2q(gate):
+        if cz_dec.count_2q(gate, backend) < iswap_dec.count_2q(gate, backend):
             return cz_dec(gate)
-        if cz_dec.count_2q(gate) > iswap_dec.count_2q(gate):
-            return iswap_dec(gate)
+        if cz_dec.count_2q(gate, backend) > iswap_dec.count_2q(gate, backend):
+            return iswap_dec(gate, backend)
         # If equal check the decomposition with less 1 qubit gates.
         # This is never used for now but may be useful for future generalization
-        if cz_dec.count_1q(gate) < iswap_dec.count_1q(gate):  # pragma: no cover
-            return cz_dec(gate)
-        return iswap_dec(gate)  # pragma: no cover
+        if cz_dec.count_1q(gate, backend) < iswap_dec.count_1q(
+            gate, backend
+        ):  # pragma: no cover
+            return cz_dec(gate, backend)
+        return iswap_dec(gate, backend)  # pragma: no cover
 
     if native_gates & NativeGates.CZ:
-        return cz_dec(gate)
+        return cz_dec(gate, backend)
 
     if native_gates & NativeGates.iSWAP:
         if gate.__class__ in iswap_dec.decompositions:
-            return iswap_dec(gate)
+            return iswap_dec(gate, backend)
 
         # First decompose into CZ
-        cz_decomposed = cz_dec(gate)
+        cz_decomposed = cz_dec(gate, backend)
         # Then CZ are decomposed into iSWAP
         iswap_decomposed = []
         for g in cz_decomposed:
             # Need recursive function as gates.Unitary is not in iswap_dec
-            for g_translated in translate_gate(g, native_gates=native_gates):
+            for g_translated in translate_gate(
+                g, native_gates=native_gates, backend=backend
+            ):
                 iswap_decomposed.append(g_translated)
         return iswap_decomposed
 
@@ -253,7 +270,7 @@ def _translate_two_qubit_gates(gate: gates.Gate, native_gates: NativeGates):
     # No CZ, iSWAP gates in the native gate set
     # Decompose CNOT, CZ, SWAP gates into CNOT gates
     if native_gates & NativeGates.CNOT:
-        return cnot_dec_temp(gate)
+        return cnot_dec_temp(gate, backend)
 
     raise_error(
         DecompositionError,
diff --git a/tests/test_backends.py b/tests/test_backends.py
index 1be00e6d22..0c855c9686 100644
--- a/tests/test_backends.py
+++ b/tests/test_backends.py
@@ -130,19 +130,3 @@ def test_list_available_backends():
     assert available_backends == list_available_backends(
         "qibojit", "qibolab", "qibo-cloud-backends", "qibotn"
     )
-
-
-def test_gradients_pytorch():
-    from qibo.backends import PyTorchBackend  # pylint: disable=import-outside-toplevel
-
-    backend = PyTorchBackend()
-    gate = gates.RX(0, 0.1)
-    matrix = gate.matrix(backend)
-    assert matrix.requires_grad
-    assert backend.gradients
-    backend.requires_grad(False)
-    gate = gates.RX(0, 0.1)
-    matrix = gate.matrix(backend)
-    assert not matrix.requires_grad
-    assert not backend.gradients
-    assert not backend.matrices.requires_grad
diff --git a/tests/test_backends_torch_gradients.py b/tests/test_backends_torch_gradients.py
new file mode 100644
index 0000000000..08bacb44af
--- /dev/null
+++ b/tests/test_backends_torch_gradients.py
@@ -0,0 +1,91 @@
+import sys
+
+import numpy as np
+import pytest
+
+from qibo import gates, models
+from qibo.backends import PyTorchBackend
+from qibo.quantum_info import infidelity
+
+
+def test_torch_gradients():
+    backend = PyTorchBackend()
+    backend.np.manual_seed(42)
+    nepochs = 400
+    optimizer = backend.np.optim.Adam
+    target_state = backend.np.rand(4, dtype=backend.np.complex128)
+    target_state = target_state / backend.np.norm(target_state)
+    params = backend.np.rand(4, dtype=backend.np.float64, requires_grad=True)
+    c = models.Circuit(2)
+    c.add(gates.RX(0, params[0]))
+    c.add(gates.RY(1, params[1]))
+    c.add(gates.U2(1, params[2], params[3]))
+
+    initial_params = params.clone()
+    initial_loss = infidelity(
+        target_state, backend.execute_circuit(c).state(), backend=backend
+    )
+
+    optimizer = optimizer([params], lr=0.01)
+    for _ in range(nepochs):
+        optimizer.zero_grad()
+        c.set_parameters(params)
+        final_state = backend.execute_circuit(c).state()
+        loss = infidelity(target_state, final_state, backend=backend)
+        loss.backward()
+        grad = params.grad.clone().norm()
+        optimizer.step()
+
+    assert initial_loss > loss
+    assert initial_params[0] != params[0]
+    assert grad.item() < 10e-3
+
+
+@pytest.mark.skipif(
+    sys.platform != "linux", reason="Tensorflow available only when testing on linux."
+)
+def test_torch_tensorflow_gradients():
+    backend = PyTorchBackend()
+
+    import tensorflow as tf  # pylint: disable=import-outside-toplevel
+
+    from qibo.backends.tensorflow import (  # pylint: disable=import-outside-toplevel
+        TensorflowBackend,
+    )
+
+    target_state = backend.np.tensor([0.0, 1.0], dtype=backend.np.complex128)
+    param = backend.np.tensor([0.1], dtype=backend.np.float64, requires_grad=True)
+    c = models.Circuit(1)
+    c.add(gates.RX(0, param[0]))
+
+    optimizer = backend.np.optim.SGD
+    optimizer = optimizer([param], lr=1)
+    c.set_parameters(param)
+    final_state = backend.execute_circuit(c).state()
+    loss = infidelity(target_state, final_state, backend=backend)
+    loss.backward()
+    torch_param_grad = param.grad.clone().item()
+    optimizer.step()
+    torch_param = param.clone().item()
+
+    tf_backend = TensorflowBackend()
+
+    target_state = tf.constant([0.0, 1.0], dtype=tf.complex128)
+    param = tf.Variable([0.1], dtype=tf.float64)
+    c = models.Circuit(1)
+    c.add(gates.RX(0, param[0]))
+
+    optimizer = tf.optimizers.SGD(learning_rate=1.0)
+
+    with tf.GradientTape() as tape:
+        c.set_parameters(param)
+        final_state = tf_backend.execute_circuit(c).state()
+        loss = infidelity(target_state, final_state, backend=tf_backend)
+
+    grads = tape.gradient(loss, [param])
+    tf_param_grad = grads[0].numpy()[0]
+    optimizer.apply_gradients(zip(grads, [param]))
+    tf_param = param.numpy()[0]
+
+    assert np.allclose(torch_param_grad, tf_param_grad, atol=1e-7, rtol=1e-7)
+    assert np.allclose(torch_param, tf_param, atol=1e-7, rtol=1e-7)
diff --git a/tests/test_cirq.py b/tests/test_cirq.py
index 8ff8498d7a..0a99e1405e 100644
--- a/tests/test_cirq.py
+++ b/tests/test_cirq.py
@@ -60,9 +60,13 @@ def assert_gates_equivalent(
     assert c.depth == target_depth
     if accelerators and not backend.supports_multigpu:
         with pytest.raises(NotImplementedError):
-            final_state = backend.execute_circuit(c, np.copy(initial_state)).state()
+            final_state = backend.execute_circuit(
+                c, backend.cast(initial_state, copy=True)
+            ).state()
     else:
-        final_state = backend.execute_circuit(c, np.copy(initial_state)).state()
+        final_state = backend.execute_circuit(
+            c, backend.cast(initial_state, copy=True)
+        ).state()
         backend.assert_allclose(final_state, target_state, atol=atol)
 
 
diff --git a/tests/test_gates_gates.py b/tests/test_gates_gates.py
index 153528cd17..f47f26d613 100644
--- a/tests/test_gates_gates.py
+++ b/tests/test_gates_gates.py
@@ -1382,8 +1382,8 @@ def test_unitary_initialization(backend):
 def test_unitary_common_gates(backend):
     target_state = apply_gates(backend, [gates.X(0), gates.H(1)], nqubits=2)
     gatelist = [
-        gates.Unitary(np.array([[0, 1], [1, 0]]), 0),
-        gates.Unitary(np.array([[1, 1], [1, -1]]) / np.sqrt(2), 1),
+        gates.Unitary(backend.cast([[0.0, 1.0], [1.0, 0.0]]), 0),
+        gates.Unitary(backend.cast([[1.0, 1.0], [1.0, -1.0]]) / np.sqrt(2), 1),
     ]
     final_state = apply_gates(backend, gatelist, nqubits=2)
     backend.assert_allclose(final_state, target_state, atol=1e-6)
@@ -1405,7 +1405,7 @@ def test_unitary_common_gates(backend):
             [np.sin(thetay / 2), np.cos(thetay / 2)],
         ]
     )
-    cnot = np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]])
+    cnot = np.array([[1.0, 0, 0, 0], [0, 1.0, 0, 0], [0, 0, 0, 1.0], [0, 0, 1.0, 0]])
     gatelist = [gates.Unitary(rx, 0), gates.Unitary(ry, 1), gates.Unitary(cnot, 0, 1)]
     final_state = apply_gates(backend, gatelist, nqubits=2)
     backend.assert_allclose(final_state, target_state, atol=1e-6)
diff --git a/tests/test_models_circuit.py b/tests/test_models_circuit.py
index 8ffa586740..c6a72a39ff 100644
--- a/tests/test_models_circuit.py
+++ b/tests/test_models_circuit.py
@@ -351,19 +351,19 @@ def test_circuit_light_cone():
 OPENQASM 2.0;
 include "qelib1.inc";
 qreg q[6];
-ry(0) q[0];
-ry(0) q[1];
-ry(0) q[2];
-ry(0) q[3];
-ry(0) q[4];
-ry(0) q[5];
+ry(0.0) q[0];
+ry(0.0) q[1];
+ry(0.0) q[2];
+ry(0.0) q[3];
+ry(0.0) q[4];
+ry(0.0) q[5];
 cz q[0],q[1];
 cz q[2],q[3];
 cz q[4],q[5];
-ry(0) q[1];
-ry(0) q[2];
-ry(0) q[3];
-ry(0) q[4];
+ry(0.0) q[1];
+ry(0.0) q[2];
+ry(0.0) q[3];
+ry(0.0) q[4];
 cz q[1],q[2];
 cz q[3],q[4];"""
     assert qubit_map == {2: 0, 3: 1, 4: 2, 5: 3, 6: 4, 7: 5}
diff --git a/tests/test_models_circuit_parametrized.py b/tests/test_models_circuit_parametrized.py
index 1a761f16eb..4e04ccd248 100644
--- a/tests/test_models_circuit_parametrized.py
+++ b/tests/test_models_circuit_parametrized.py
@@ -39,12 +39,12 @@ def test_set_parameters_with_list(backend, trainable):
     params = [0.123, 0.456, (0.789, 0.321)]
     c = Circuit(3)
     if trainable:
-        c.add(gates.RX(0, theta=0, trainable=trainable))
+        c.add(gates.RX(0, theta=0.0, trainable=trainable))
     else:
         c.add(gates.RX(0, theta=params[0], trainable=trainable))
-    c.add(gates.RY(1, theta=0))
+    c.add(gates.RY(1, theta=0.0))
     c.add(gates.CZ(1, 2))
-    c.add(gates.fSim(0, 2, theta=0, phi=0))
+    c.add(gates.fSim(0, 2, theta=0.0, phi=0.0))
     c.add(gates.H(2))
     # execute once
     final_state = backend.execute_circuit(c)
@@ -78,16 +78,16 @@ def test_circuit_set_parameters_ungates(backend, trainable, accelerators):
         trainable_params = [0.1, 0.3, (0.4, 0.5)]
 
     c = Circuit(3, accelerators)
-    c.add(gates.RX(0, theta=0))
+    c.add(gates.RX(0, theta=0.0))
     if trainable:
-        c.add(gates.CRY(0, 1, theta=0, trainable=trainable))
+        c.add(gates.CRY(0, 1, theta=0.0, trainable=trainable))
     else:
         c.add(gates.CRY(0, 1, theta=params[1], trainable=trainable))
     c.add(gates.CZ(1, 2))
-    c.add(gates.U1(2, theta=0))
-    c.add(gates.CU2(0, 2, phi=0, lam=0))
+    c.add(gates.U1(2, theta=0.0))
+    c.add(gates.CU2(0, 2, phi=0.0, lam=0.0))
     if trainable:
-        c.add(gates.U3(1, theta=0, phi=0, lam=0, trainable=trainable))
+        c.add(gates.U3(1, theta=0.0, phi=0.0, lam=0.0, trainable=trainable))
     else:
         c.add(gates.U3(1, *params[4], trainable=trainable))
     # execute once
@@ -127,7 +127,7 @@ def test_circuit_set_parameters_with_unitary(backend, trainable, accelerators):
     """Check updating parameters of circuit that contains ``Unitary`` gate."""
     params = [0.1234, np.random.random((4, 4))]
     c = Circuit(4, accelerators)
-    c.add(gates.RX(0, theta=0))
+    c.add(gates.RX(0, theta=0.0))
     if trainable:
         c.add(gates.Unitary(np.zeros((4, 4)), 1, 2, trainable=trainable))
         trainable_params = list(params)
diff --git a/tests/test_models_encodings.py b/tests/test_models_encodings.py
index 1a935d75e3..c8f1dacf14 100644
--- a/tests/test_models_encodings.py
+++ b/tests/test_models_encodings.py
@@ -135,7 +135,11 @@ def test_unary_encoder(backend, nqubits, architecture, kind):
     indexes = np.flatnonzero(backend.to_numpy(state))
     state = backend.np.real(state[indexes])
 
-    backend.assert_allclose(state, backend.cast(data) / backend.calculate_norm(data, 2))
+    backend.assert_allclose(
+        state,
+        backend.cast(data, dtype=np.float64) / backend.calculate_norm(data, 2),
+        rtol=1e-5,
+    )
 
 
 @pytest.mark.parametrize("seed", [None, 10, np.random.default_rng(10)])
diff --git a/tests/test_models_hep.py b/tests/test_models_hep.py
index fdf0a46112..19dd64146e 100644
--- a/tests/test_models_hep.py
+++ b/tests/test_models_hep.py
@@ -108,4 +108,6 @@ def test_qpdf(backend, ansatz, layers, nqubits, multi_output, output):
     np.random.seed(0)
     params = np.random.rand(model.nparams)
     result = model.predict(params, [0.1])
-    np.testing.assert_allclose(result, output, atol=1e-5)
+    atol = 1e-5
+    rtol = 1e-7
+    np.testing.assert_allclose(result, output, rtol=rtol, atol=atol)
diff --git a/tests/test_models_qft.py b/tests/test_models_qft.py
index 90daa4fe9c..1ead4cb865 100644
--- a/tests/test_models_qft.py
+++ b/tests/test_models_qft.py
@@ -48,10 +48,10 @@ def test_qft_matrix(backend, nqubits):
     c = models.QFT(nqubits)
     dim = 2**nqubits
     target_matrix = qft_matrix(dim)
-    backend.assert_allclose(c.unitary(backend), target_matrix)
+    backend.assert_allclose(c.unitary(backend), target_matrix, atol=1e-6, rtol=1e-6)
     c = c.invert()
     target_matrix = qft_matrix(dim, inverse=True)
-    backend.assert_allclose(c.unitary(backend), target_matrix)
+    backend.assert_allclose(c.unitary(backend), target_matrix, atol=1e-6, rtol=1e-6)
 
 
 @pytest.mark.parametrize("density_matrix", [False, True])
@@ -69,7 +69,7 @@ def test_qft_execution(backend, nqubits, random, density_matrix):
 
     final_state = backend.execute_circuit(c, backend.np.copy(initial_state))._state
     target_state = exact_qft(initial_state, density_matrix, backend)
-    backend.assert_allclose(final_state, target_state)
+    backend.assert_allclose(final_state, target_state, atol=1e-6, rtol=1e-6)
 
 
 def test_qft_errors():
diff --git a/tests/test_models_variational.py b/tests/test_models_variational.py
index c596990488..5ec737610f 100644
--- a/tests/test_models_variational.py
+++ b/tests/test_models_variational.py
@@ -124,10 +124,11 @@ def test_vqe(backend, method, options, compile, filename):
         circuit.add(gates.RY(q, theta=1.0))
     hamiltonian = hamiltonians.XXZ(nqubits=nqubits, backend=backend)
     np.random.seed(0)
-    initial_parameters = np.random.uniform(0, 2 * np.pi, 2 * nqubits * layers + nqubits)
     initial_parameters = backend.cast(
-        initial_parameters, dtype=initial_parameters.dtype
+        np.random.uniform(0, 2 * np.pi, 2 * nqubits * layers + nqubits), dtype="float64"
     )
+    if backend.name == "pytorch":
+        initial_parameters.requires_grad = True
     v = models.VQE(circuit, hamiltonian)
 
     loss_values = []
@@ -150,7 +151,7 @@ def callback(parameters, loss_values=loss_values, vqe=v):
         shutil.rmtree("outcmaes")
     if filename is not None:
         assert_regression_fixture(backend, params, filename)
-    assert best == min(loss_values)
+    backend.assert_allclose(best, min(loss_values), rtol=1e-6, atol=1e-6)
 
     # test energy fluctuation
     state = backend.np.ones(2**nqubits) / np.sqrt(2**nqubits)
@@ -275,8 +276,9 @@ def test_qaoa_optimization(backend, method, options, dense, filename):
         pytest.skip("Skipping scipy optimizers for tensorflow and pytorch.")
     h = hamiltonians.XXZ(3, dense=dense, backend=backend)
     qaoa = models.QAOA(h)
-    initial_p = [0.05, 0.06, 0.07, 0.08]
-    initial_p = backend.cast(initial_p, dtype=np.float64)
+    initial_p = backend.cast([0.05, 0.06, 0.07, 0.08], dtype="float64")
+    if backend.name == "pytorch":
+        initial_p.requires_grad = True
     best, params, _ = qaoa.minimize(initial_p, method=method, options=options)
     if filename is not None:
         assert_regression_fixture(backend, params, filename)
diff --git a/tests/test_noise.py b/tests/test_noise.py
index c7518a3a40..fb634fa1cd 100644
--- a/tests/test_noise.py
+++ b/tests/test_noise.py
@@ -85,8 +85,8 @@ def test_kraus_error(backend, density_matrix, nshots):
 @pytest.mark.parametrize("density_matrix", [False, True])
 @pytest.mark.parametrize("nshots", [10, 100])
 def test_unitary_error(backend, density_matrix, nshots):
-    u1 = np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]])
-    u2 = np.array([[0, 1, 0, 0], [1, 0, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]])
+    u1 = np.array([[1.0, 0, 0, 0], [0, 1.0, 0, 0], [0, 0, 0, 1.0], [0, 0, 1.0, 0]])
+    u2 = np.array([[0, 1.0, 0, 0], [1.0, 0, 0, 0], [0, 0, 0, 1.0], [0, 0, 1.0, 0]])
     qubits = (0, 1)
     p1, p2 = (0.3, 0.7)
 
diff --git a/tests/test_quantum_info_entanglement.py b/tests/test_quantum_info_entanglement.py
index a461c16495..ed7e98f300 100644
--- a/tests/test_quantum_info_entanglement.py
+++ b/tests/test_quantum_info_entanglement.py
@@ -149,7 +149,7 @@ def test_entanglement_fidelity(backend, channel, nqubits, check_hermitian):
 
 def test_meyer_wallach_entanglement(backend):
     with pytest.raises(TypeError):
-        state = np.random.rand(2, 3, 2)
+        state = np.random.rand(2, 3, 2).astype(complex)
         state = backend.cast(state, dtype=state.dtype)
         test = meyer_wallach_entanglement(state, backend=backend)
 
@@ -165,11 +165,11 @@ def test_meyer_wallach_entanglement(backend):
     state2 = backend.execute_circuit(circuit2).state()
 
     backend.assert_allclose(
-        meyer_wallach_entanglement(state1, backend=backend), 0.0, atol=PRECISION_TOL
+        meyer_wallach_entanglement(state1, backend=backend), 0.0, atol=1e-6, rtol=1e-6
     )
 
     backend.assert_allclose(
-        meyer_wallach_entanglement(state2, backend=backend), 1, atol=PRECISION_TOL
+        meyer_wallach_entanglement(state2, backend=backend), 1.0, atol=1e-6, rtol=1e-6
     )
 
 
diff --git a/tests/test_tomography_gate_set_tomography.py b/tests/test_tomography_gate_set_tomography.py
index bbed8cff06..d6d2ffd2df 100644
--- a/tests/test_tomography_gate_set_tomography.py
+++ b/tests/test_tomography_gate_set_tomography.py
@@ -285,7 +285,7 @@ def test_GST_with_transpiler(backend):
             Preprocessing(connectivity),
             Random(connectivity),
             Sabre(connectivity),
-            Unroller(NativeGates.default()),
+            Unroller(NativeGates.default(), backend=backend),
         ],
         int_qubit_names=True,
     )
diff --git a/tests/test_transpiler_decompositions.py b/tests/test_transpiler_decompositions.py
index 663f73f5ff..44c299cc3b 100644
--- a/tests/test_transpiler_decompositions.py
+++ b/tests/test_transpiler_decompositions.py
@@ -21,7 +21,7 @@ def assert_matrices_allclose(gate, natives, backend):
     target_unitary = target_matrix / normalisation
 
     circuit = Circuit(len(gate.qubits))
-    circuit.add(translate_gate(gate, natives))
+    circuit.add(translate_gate(gate, natives, backend=backend))
     native_matrix = circuit.unitary(backend)
     # Remove global phase from native matrix
     normalisation = np.power(
@@ -34,10 +34,10 @@ def assert_matrices_allclose(gate, natives, backend):
     # There can still be phase differences of -1, -1j, 1j
     c = 0
     for phase in [1, -1, 1j, -1j]:
-        if np.allclose(
-            backend.to_numpy(phase * native_unitary),
-            backend.to_numpy(target_unitary),
-            atol=1e-12,
+        if backend.np.allclose(
+            phase * native_unitary,
+            target_unitary,
+            atol=1e-6,
         ):
             c = 1
     backend.assert_allclose(c, 1)
@@ -229,15 +229,15 @@ def test_unitary_to_native(backend, nqubits, natives_1q, natives_2q, seed):
     assert_matrices_allclose(gate, natives_1q | natives_2q | default_natives, backend)
 
 
-def test_count_1q():
+def test_count_1q(backend):
     from qibo.transpiler.unroller import cz_dec
 
-    np.testing.assert_allclose(cz_dec.count_1q(gates.CNOT(0, 1)), 2)
-    np.testing.assert_allclose(cz_dec.count_1q(gates.CRX(0, 1, 0.1)), 2)
+    np.testing.assert_allclose(cz_dec.count_1q(gates.CNOT(0, 1), backend), 2)
+    np.testing.assert_allclose(cz_dec.count_1q(gates.CRX(0, 1, 0.1), backend), 2)
 
 
-def test_count_2q():
+def test_count_2q(backend):
     from qibo.transpiler.unroller import cz_dec
 
-    np.testing.assert_allclose(cz_dec.count_2q(gates.CNOT(0, 1)), 1)
-    np.testing.assert_allclose(cz_dec.count_2q(gates.CRX(0, 1, 0.1)), 2)
+    np.testing.assert_allclose(cz_dec.count_2q(gates.CNOT(0, 1), backend), 1)
+    np.testing.assert_allclose(cz_dec.count_2q(gates.CRX(0, 1, 0.1), backend), 2)
diff --git a/tests/test_transpiler_unitary_decompositions.py b/tests/test_transpiler_unitary_decompositions.py
index df641208bb..49f8cf2fe9 100644
--- a/tests/test_transpiler_unitary_decompositions.py
+++ b/tests/test_transpiler_unitary_decompositions.py
@@ -123,8 +123,8 @@ def test_ud_eigenvalues(backend, seed):
             @ backend.cast(bell_basis)
         )
         ud_diag = backend.np.diag(ud_bell)
-        backend.assert_allclose(backend.np.diag(ud_diag), ud_bell, atol=PRECISION_TOL)
-        backend.assert_allclose(backend.np.prod(ud_diag), 1)
+        backend.assert_allclose(backend.np.diag(ud_diag), ud_bell, atol=1e-6, rtol=1e-6)
+        backend.assert_allclose(backend.np.prod(ud_diag), 1, atol=1e-6, rtol=1e-6)
 
 
 @pytest.mark.parametrize("seed", [None, 10, np.random.default_rng(10)])
@@ -139,7 +139,7 @@ def test_calculate_h_vector(backend, seed):
         assert ud_diag is not None
         hx, hy, hz = calculate_h_vector(ud_diag, backend=backend)
         target_matrix = bell_unitary(hx, hy, hz, backend)
-        backend.assert_allclose(ud, target_matrix, atol=PRECISION_TOL)
+        backend.assert_allclose(ud, target_matrix, atol=1e-6, rtol=1e-6)
 
 
 def test_cnot_decomposition(backend):
@@ -148,7 +148,7 @@ def test_cnot_decomposition(backend):
     c = Circuit(2)
     c.add(cnot_decomposition(0, 1, hx, hy, hz, backend))
     final_matrix = c.unitary(backend)
-    backend.assert_allclose(final_matrix, target_matrix, atol=PRECISION_TOL)
+    backend.assert_allclose(final_matrix, target_matrix, atol=1e-6, rtol=1e-6)
 
 
 def test_cnot_decomposition_light(backend):
@@ -157,7 +157,7 @@ def test_cnot_decomposition_light(backend):
     c = Circuit(2)
     c.add(cnot_decomposition_light(0, 1, hx, hy, backend))
     final_matrix = c.unitary(backend)
-    backend.assert_allclose(final_matrix, target_matrix, atol=PRECISION_TOL)
+    backend.assert_allclose(final_matrix, target_matrix, atol=1e-6, rtol=1e-6)
 
 
 @pytest.mark.parametrize("seed", [None, 10, np.random.default_rng(10)])
@@ -170,7 +170,7 @@ def test_two_qubit_decomposition(backend, seed):
         c = Circuit(2)
         c.add(two_qubit_decomposition(0, 1, unitary, backend=backend))
         final_matrix = c.unitary(backend)
-        backend.assert_allclose(final_matrix, unitary, atol=PRECISION_TOL)
+        backend.assert_allclose(final_matrix, unitary, atol=1e-6, rtol=1e-6)
 
 
 @pytest.mark.parametrize("gatename", ["CNOT", "CZ", "SWAP", "iSWAP", "fSim", "I"])
@@ -191,7 +191,7 @@ def test_two_qubit_decomposition_common_gates(backend, gatename):
         c = Circuit(2)
         c.add(two_qubit_decomposition(0, 1, matrix, backend=backend))
         final_matrix = c.unitary(backend)
-        backend.assert_allclose(final_matrix, matrix, atol=PRECISION_TOL)
+        backend.assert_allclose(final_matrix, matrix, atol=1e-6, rtol=1e-6)
 
 
 @pytest.mark.parametrize("hz_zero", [False, True])
@@ -203,7 +203,7 @@ def test_two_qubit_decomposition_bell_unitary(backend, hz_zero):
     c = Circuit(2)
     c.add(two_qubit_decomposition(0, 1, unitary, backend=backend))
     final_matrix = c.unitary(backend)
-    backend.assert_allclose(final_matrix, unitary, atol=PRECISION_TOL)
+    backend.assert_allclose(final_matrix, unitary, atol=1e-6, rtol=1e-6)
 
 
 def test_two_qubit_decomposition_no_entanglement(backend):
@@ -224,4 +224,4 @@ def test_two_qubit_decomposition_no_entanglement(backend):
     else:
         c.add(two_qubit_decomposition(0, 1, matrix, backend=backend))
         final_matrix = c.unitary(backend)
-        backend.assert_allclose(final_matrix, matrix, atol=PRECISION_TOL)
+        backend.assert_allclose(final_matrix, matrix, atol=1e-6, rtol=1e-6)