multi-qubit Pauli noise channel

doc/source/api-reference/qibo.rst

@@ -614,6 +614,13 @@ Pauli noise channel
     :members:
     :member-order: bysource
 
+Generalized Pauli noise channel
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+.. autoclass:: qibo.gates.GeneralizedPauliNoiseChannel
+    :members:
+    :member-order: bysource
+
 Depolarizing channel
 ^^^^^^^^^^^^^^^^^^^^
 

src/qibo/gates/channels.py

@@ -1,4 +1,6 @@
+import warnings
 from itertools import product
+from typing import Tuple
 
 from qibo.config import PRECISION_TOL, raise_error
 from qibo.gates.abstract import Gate
@@ -88,9 +90,8 @@ def __init__(self, ops):
                 if shape != (rank, rank):
                     raise_error(
                         ValueError,
-                        "Invalid Krauss operator shape {} for "
-                        "acting on {} qubits."
-                        "".format(shape, len(qubits)),
+                        f"Invalid Krauss operator shape {shape} for "
+                        + f"acting on {len(qubits)} qubits.",
                     )
                 qubitset.update(qubits)
                 gates.append(Unitary(matrix, *list(qubits)))
@@ -243,9 +244,8 @@ def __init__(self, probabilities, ops):
         if len(probabilities) != len(ops):
             raise_error(
                 ValueError,
-                "Probabilities list has length {} while "
-                "{} gates were given."
-                "".format(len(probabilities), len(ops)),
+                f"Probabilities list has length {len(probabilities)} while "
+                + f"{len(ops)} gates were given.",
             )
         for p in probabilities:
             if p < 0 or p > 1:
@@ -292,6 +292,12 @@ class PauliNoiseChannel(UnitaryChannel):
     """
 
     def __init__(self, q, px=0, py=0, pz=0):
+        warnings.warn(
+            "This channel will be removed in a later release. "
+            + "Use GeneralizedPauliNoiseChannel instead.",
+            DeprecationWarning,
+        )
+
         probs, gates = [], []
         for p, gate in [(px, X), (py, Y), (pz, Z)]:
             if p > 0:
@@ -306,13 +312,92 @@ def __init__(self, q, px=0, py=0, pz=0):
         self.init_kwargs = {"px": px, "py": py, "pz": pz}
 
 
+class GeneralizedPauliNoiseChannel(UnitaryChannel):
+    """Multi-qubit noise channel that applies Pauli operators with given probabilities.
+
+    Implements the following transformation:
+
+    .. math::
+        \\mathcal{E}(\\rho ) = \\left (1 - \\sum _{k} p_{k} \\right ) \\, \\rho +
+                                \\sum_{k} \\, p_{k} \\, P_{k} \\, \\rho \\, P_{k}
+
+
+    where :math:`P_{k}` is the :math:`k`-th Pauli ``string`` and :math:`p_{k}` is
+    the probability associated to :math:`P_{k}`.
+
+    Example:
+        .. testcode::
+
+            import numpy as np
+
+            from itertools import product
+
+            from qibo.gates.channels import GeneralizedPauliNoiseChannel
+
+            qubits = (0, 2)
+            nqubits = len(qubits)
+
+            # excluding the Identity operator
+            paulis = list(product(["I", "X"], repeat=nqubits))[1:]
+            # this next line is optional
+            paulis = [''.join(pauli) for pauli in paulis]
+
+            probabilities = np.random.rand(len(paulis) + 1)
+            probabilities /= np.sum(probabilities)
+            #Excluding probability of Identity operator
+            probabilities = probabilities[1:]
+
+            channel = GeneralizedPauliNoiseChannel(
+                qubits, list(zip(paulis, probabilities))
+            )
+
+    This channel can be simulated using either density matrices or state vectors
+    and sampling with repeated execution.
+    See :ref:`How to perform noisy simulation? <noisy-example>` for more
+    information.
+
+    Args:
+        qubits (int or list or tuple): Qubits that the noise acts on.
+        operators (list): list of operators as pairs :math:`(P_{k}, p_{k})`.
+    """
+
+    def __init__(self, qubits: Tuple[int, list, tuple], operators: list):
+        warnings.warn(
+            "The class GeneralizedPauliNoiseChannel will be renamed "
+            + "PauliNoiseChannel in a later release."
+        )
+
+        if isinstance(qubits, int) is True:
+            qubits = (qubits,)
+
+        probabilities, paulis = [], []
+        for pauli, probability in operators:
+            probabilities.append(probability)
+            paulis.append(pauli)
+
+        single_paulis = {"I": I, "X": X, "Y": Y, "Z": Z}
+
+        gates = []
+        for pauli in paulis:
+            fgate = FusedGate(*qubits)
+            for q, p in zip(qubits, pauli):
+                fgate.append(single_paulis[p](q))
+            gates.append(fgate)
+        self.gates = tuple(gates)
+        self.coefficients = tuple(probabilities)
+
+        super().__init__(probabilities, gates)
+        self.name = "GeneralizedPauliNoiseChannel"
+
+
 class DepolarizingChannel(Channel):
     """:math:`n`-qubit Depolarizing quantum error channel,
 
     .. math::
         \\mathcal{E}(\\rho ) = (1 - \\lambda) \\rho +\\lambda \\text{Tr}_q[\\rho]\\otimes \\frac{I}{2^n}
 
-    where :math:`\\lambda` is the depolarizing error parameter and :math:`0 \\le \\lambda \\le 4^n / (4^n - 1)`.
+    where :math:`\\lambda` is the depolarizing error parameter
+    and :math:`0 \\le \\lambda \\le 4^n / (4^n - 1)`.
 
     * If :math:`\\lambda = 1` this is a completely depolarizing channel
       :math:`E(\\rho) = I / 2^n`
@@ -325,7 +410,7 @@ class DepolarizingChannel(Channel):
         lam (float): Depolarizing error parameter.
     """
 
-    def __init__(self, q, lam=0):
+    def __init__(self, q, lam: str = 0):
         super().__init__()
         num_qubits = len(q)
         num_terms = 4**num_qubits
@@ -344,9 +429,11 @@ def __init__(self, q, lam=0):
 
     def apply_density_matrix(self, backend, state, nqubits):
         lam = self.init_kwargs["lam"]
-        return (1 - lam) * backend.cast(state) + lam / 2**nqubits * backend.cast(
-            I(*range(nqubits)).asmatrix(backend)
-        )
+        state_evolved = (1 - lam) * backend.cast(state) + (
+            lam / 2**nqubits
+        ) * backend.cast(I(*range(nqubits)).asmatrix(backend))
+
+        return state_evolved
 
     def apply(self, backend, state, nqubits):
         num_qubits = len(self.target_qubits)
@@ -355,12 +442,13 @@ def apply(self, backend, state, nqubits):
         probs = (num_terms - 1) * [prob_pauli]
         gates = []
         for pauli_list in list(product([I, X, Y, Z], repeat=num_qubits))[1::]:
-            fgate = FusedGate(*range(num_qubits))
+            fgate = FusedGate(*self.target_qubits)
             for j, pauli in enumerate(pauli_list):
                 fgate.append(pauli(j))
-            gates.append(Unitary(backend.asmatrix_fused(fgate), *self.target_qubits))
+            gates.append(fgate)
         self.gates = tuple(gates)
         self.coefficients = tuple(probs)
+
         return backend.apply_channel(self, state, nqubits)
 
 

src/qibo/quantum_info/metrics.py

@@ -328,7 +328,6 @@ def gate_error(channel, target=None):
 
     Returns:
         float: Gate error between :math:`\\mathcal{E}` and :math:`\\mathcal{U}`.
-
     """
 
     return 1 - average_gate_fidelity(channel, target)

src/qibo/tests/test_gates_channels.py

@@ -159,6 +159,32 @@ def test_pauli_noise_channel(backend):
     assert norm < PRECISION_TOL
 
 
+def test_generalized_pauli_noise_channel(backend):
+    initial_rho = random_density_matrix(2**2)
+    qubits = (1,)
+    channel = gates.GeneralizedPauliNoiseChannel(qubits, [("X", 0.3)])
+    final_rho = backend.apply_channel_density_matrix(channel, np.copy(initial_rho), 2)
+    gate = gates.X(1)
+    target_rho = backend.apply_gate_density_matrix(gate, np.copy(initial_rho), 2)
+    target_rho = 0.3 * backend.to_numpy(target_rho)
+    target_rho += 0.7 * initial_rho
+    backend.assert_allclose(final_rho, target_rho)
+
+    basis = ["X", "Y", "Z"]
+    pnp = np.array([0.1, 0.02, 0.05])
+    a0 = 1
+    a1 = 1 - 2 * pnp[1] - 2 * pnp[2]
+    a2 = 1 - 2 * pnp[0] - 2 * pnp[2]
+    a3 = 1 - 2 * pnp[0] - 2 * pnp[1]
+    test_representation = np.diag([a0, a1, a2, a3])
+
+    liouville = gates.GeneralizedPauliNoiseChannel(
+        0, list(zip(basis, pnp))
+    ).to_pauli_liouville(True, backend)
+    norm = np.linalg.norm(backend.to_numpy(liouville) - test_representation)
+    assert norm < PRECISION_TOL
+
+
 def test_depolarizing_channel(backend):
     initial_rho = random_density_matrix(2**3)
     lam = 0.3