Qutrit channel depolarizing

doc/introduction/operations.rst

@@ -495,6 +495,21 @@ Qutrit State preparation
 
 :html:`</div>`
 
+.. _intro_ref_ops_qutrit_channels:
+
+Qutrit noisy channels
+^^^^^^^^^^^^^^^^^^^^^^^^
+
+
+:html:`<div class="summary-table">`
+
+.. autosummary::
+    :nosignatures:
+
+    ~pennylane.QutritDepolarizingChannel
+
+:html:`</div>`
+
 .. _intro_ref_ops_qutrit_obs:
 
 Qutrit Observables

doc/releases/changelog-dev.md

@@ -98,6 +98,11 @@
 * ``qml.load`` has been removed in favour of more specific functions, such as ``qml.from_qiskit``, etc.
   [(#5654)](https://github.com/PennyLaneAI/pennylane/pull/5654)
 
+<h4>Community contributions 🥳</h4>
+
+* ``qml.QutritDepolarizingChannel`` has been added, allowing for depolarizing noise to be simulated on the `default.qutrit.mixed` device.
+  [(#5502)](https://github.com/PennyLaneAI/pennylane/pull/5502)
+
 <h3>Deprecations 👋</h3>
 
 * ``qml.transforms.map_batch_transform`` is deprecated, since a transform can be applied directly to a batch of tapes.
@@ -131,6 +136,7 @@
 
 This release contains contributions from (in alphabetical order):
 
+Gabriel Bottrill,
 Isaac De Vlugt,
 Pietropaolo Frisoni,
 Soran Jahangiri,

pennylane/devices/default_qutrit_mixed.py

@@ -13,7 +13,6 @@
 # limitations under the License.
 """The default.qutrit.mixed device is PennyLane's standard qutrit simulator for mixed-state
 computations."""
-
 import inspect
 import logging
 from dataclasses import replace
@@ -22,6 +21,7 @@
 import numpy as np
 
 import pennylane as qml
+from pennylane.ops import _qutrit__channel__ops__ as channels
 from pennylane.tape import QuantumTape
 from pennylane.transforms.core import TransformProgram
 from pennylane.typing import Result, ResultBatch
@@ -49,7 +49,6 @@
 # always a function from a resultbatch to either a result or a result batch
 PostprocessingFn = Callable[[ResultBatch], Result_or_ResultBatch]
 
-channels = set()
 observables = {
     "THermitian",
     "GellMann",

pennylane/ops/__init__.py

@@ -37,6 +37,7 @@
 from .qutrit import __all__ as _qutrit__all__
 from .qutrit import __obs__ as _qutrit__obs__
 from .qutrit import __ops__ as _qutrit__ops__
+from .qutrit import __channels__ as _qutrit__channel__ops__
 
 _qubit__ops__ = _qubit__ops__ | _controlled_qubit__ops__
 _qubit__all__ = _qubit__all__ + list(_controlled_qubit__ops__)

pennylane/ops/qutrit/__init__.py

@@ -26,6 +26,7 @@
 from .observables import *
 from .parametric_ops import *
 from .state_preparation import *
+from .channel import QutritDepolarizingChannel
 
 # TODO: Change `qml.Identity` for qutrit support or add `qml.TIdentity` for qutrits
 __ops__ = {
@@ -47,5 +48,8 @@
     "THermitian",
     "GellMann",
 }
+__channels__ = {
+    "QutritDepolarizingChannel",
+}
 
-__all__ = list(__ops__ | __obs__)
+__all__ = list(__ops__ | __obs__ | __channels__)

pennylane/ops/qutrit/channel.py

@@ -0,0 +1,224 @@
+# Copyright 2018-2024 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+# pylint: disable=too-many-arguments
+"""
+This module contains the available built-in noisy qutrit
+quantum channels supported by PennyLane, as well as their conventions.
+"""
+import numpy as np
+
+from pennylane import math
+from pennylane.operation import Channel
+
+QUDIT_DIM = 3
+
+
+class QutritDepolarizingChannel(Channel):
+    r"""
+    Single-qutrit symmetrically depolarizing error channel.
+    This channel is modelled by the Kraus matrices generated by the following relationship:
+
+    .. math::
+        K_0 = K_{0,0} = \sqrt{1-p} \begin{bmatrix}
+                1 & 0 & 0\\
+                0 & 1 & 0\\
+                0 & 0 & 1
+                \end{bmatrix}, \quad
+        K_{i,j} = \sqrt{\frac{p}{8}}X^iZ^j
+
+    Where:
+
+    .. math::
+        X = \begin{bmatrix}
+                0 & 1 & 0 \\
+                0 & 0 & 1 \\
+                1 & 0 & 0
+                \end{bmatrix}, \quad
+        Z = \begin{bmatrix}
+                1 & 0 & 0\\
+                0 & \omega & 0\\
+                0 & 0 & \omega^2
+                \end{bmatrix}
+
+
+    These relations create the following Kraus matrices:
+
+    .. math::
+        \begin{matrix}
+            K_0 = K_{0,0} = \sqrt{1-p} \begin{bmatrix}
+                1 & 0 & 0\\
+                0 & 1 & 0\\
+                0 & 0 & 1
+                \end{bmatrix}&
+            K_1 = K_{0,1} = \sqrt{\frac{p}{8}}\begin{bmatrix}
+                1 & 0 & 0\\
+                0 & \omega & 0\\
+                0 & 0 & \omega^2
+                \end{bmatrix}&
+            K_2 = K_{0,2} = \sqrt{\frac{p}{8}}\begin{bmatrix}
+                1 & 0 & 0\\
+                0 & \omega^2 & 0\\
+                0 & 0 & \omega
+                \end{bmatrix}\\
+            K_3 = K_{1,0} = \sqrt{\frac{p}{8}}\begin{bmatrix}
+                0 & 1 & 0 \\
+                0 & 0 & 1 \\
+                1 & 0 & 0
+                \end{bmatrix}&
+            K_4 = K_{1,1} = \sqrt{\frac{p}{8}}\begin{bmatrix}
+                0 & \omega & 0 \\
+                0 & 0 & \omega^2 \\
+                1 & 0 & 0
+                \end{bmatrix}&
+            K_5 = K_{1,2} = \sqrt{\frac{p}{8}}\begin{bmatrix}
+                0 & \omega^2 & 0 \\
+                0 & 0 & \omega \\
+                1 & 0 & 0
+                \end{bmatrix}\\
+            K_6 = K_{2,0} = \sqrt{\frac{p}{8}}\begin{bmatrix}
+                0 & 0 & 1 \\
+                1 & 0 & 0 \\
+                0 & 1 & 0
+                \end{bmatrix}&
+            K_7 = K_{2,1} = \sqrt{\frac{p}{8}}\begin{bmatrix}
+                0 & 0 & \omega^2 \\
+                1 & 0 & 0 \\
+                0 & \omega & 0
+                \end{bmatrix}&
+            K_8 = K_{2,2} = \sqrt{\frac{p}{8}}\begin{bmatrix}
+                0 & 0 & \omega \\
+                1 & 0 & 0 \\
+                0 & \omega^2 & 0
+                \end{bmatrix}
+        \end{matrix}
+
+    Where :math:`\omega=\exp(\frac{2\pi}{3})`  is the third root of unity,
+    and :math:`p \in [0, 1]` is the depolarization probability, equally
+    divided in the application of all qutrit Pauli operators.
+
+    .. note::
+
+        The Kraus operators :math:`\{K_0 \ldots K_8\}` used are the representations of the single qutrit Pauli group.
+        These Pauli group operators are defined in [`1 <https://doi.org/10.48550/arXiv.quant-ph/9802007>`_] (Eq. 5).
+        The Kraus Matrices we use are adapted from [`2 <https://doi.org/10.48550/arXiv.1905.10481>`_] (Eq. 5).
+        For this definition, please make a note of the following:
+
+        * For :math:`p = 0`, the channel will be an Identity channel, i.e., a noise-free channel.
+        * For :math:`p = \frac{8}{9}`, the channel will be a fully depolarizing channel.
+        * For :math:`p = 1`, the channel will be a uniform error channel.
+
+    **Details:**
+
+    * Number of wires: 1
+    * Number of parameters: 1
+
+    Args:
+        p (float): Each qutrit Pauli operator is applied with probability :math:`\frac{p}{8}`
+        wires (Sequence[int] or int): the wire the channel acts on
+        id (str or None): String representing the operation (optional)
+    """
+
+    num_params = 1
+    num_wires = 1
+    grad_method = "A"
+    grad_recipe = ([[1, 0, 1], [-1, 0, 0]],)
+
+    def __init__(self, p, wires, id=None):
+        super().__init__(p, wires=wires, id=id)
+
+    @staticmethod
+    def compute_kraus_matrices(p):  # pylint:disable=arguments-differ
+        r"""Kraus matrices representing the qutrit depolarizing channel.
+
+         Args:
+             p (float): each qutrit Pauli gate is applied with probability :math:`\frac{p}{8}`
+
+         Returns:
+             list (array): list of Kraus matrices
+
+         **Example**
+
+         >>> np.round(qml.QutritDepolarizingChannel.compute_kraus_matrices(0.5), 3)
+         array([[[ 0.707+0.j   ,  0.   +0.j   ,  0.   +0.j   ],
+         [ 0.   +0.j   ,  0.707+0.j   ,  0.   +0.j   ],
+         [ 0.   +0.j   ,  0.   +0.j   ,  0.707+0.j   ]],
+
+        [[ 0.25 +0.j   ,  0.   +0.j   ,  0.   +0.j   ],
+         [ 0.   +0.j   , -0.125+0.217j,  0.   +0.j   ],
+         [ 0.   +0.j   ,  0.   +0.j   , -0.125-0.217j]],
+
+        [[ 0.25 +0.j   ,  0.   +0.j   ,  0.   +0.j   ],
+         [ 0.   +0.j   , -0.125-0.217j,  0.   +0.j   ],
+         [ 0.   +0.j   ,  0.   +0.j   , -0.125+0.217j]],
+
+        [[ 0.   +0.j   ,  0.25 +0.j   ,  0.   +0.j   ],
+         [ 0.   +0.j   ,  0.   +0.j   ,  0.25 +0.j   ],
+         [ 0.25 +0.j   ,  0.   +0.j   ,  0.   +0.j   ]],
+
+        [[ 0.   +0.j   , -0.125+0.217j,  0.   +0.j   ],
+         [ 0.   +0.j   ,  0.   +0.j   , -0.125-0.217j],
+         [ 0.25 +0.j   ,  0.   +0.j   ,  0.   +0.j   ]],
+
+        [[ 0.   +0.j   , -0.125-0.217j,  0.   +0.j   ],
+         [ 0.   +0.j   ,  0.   +0.j   , -0.125+0.217j],
+         [ 0.25 +0.j   ,  0.   +0.j   ,  0.   +0.j   ]],
+
+        [[ 0.   +0.j   ,  0.   +0.j   ,  0.25 +0.j   ],
+         [ 0.25 +0.j   ,  0.   +0.j   ,  0.   +0.j   ],
+         [ 0.   +0.j   ,  0.25 +0.j   ,  0.   +0.j   ]],
+
+        [[ 0.   +0.j   ,  0.   +0.j   , -0.125-0.217j],
+         [ 0.25 +0.j   ,  0.   +0.j   ,  0.   +0.j   ],
+         [ 0.   +0.j   , -0.125+0.217j,  0.   +0.j   ]],
+
+        [[ 0.   +0.j   ,  0.   +0.j   , -0.125+0.217j],
+         [ 0.25 +0.j   ,  0.   +0.j   ,  0.   +0.j   ],
+         [ 0.   +0.j   , -0.125-0.217j,  0.   +0.j   ]]])
+        """
+        if not math.is_abstract(p) and not 0.0 <= p <= 1.0:
+            raise ValueError("p must be in the interval [0,1]")
+
+        interface = math.get_interface(p)
+
+        w = math.exp(2j * np.pi / 3)
+        one = 1
+        z = 0
+
+        if interface == "tensorflow":
+            p = math.cast_like(p, 1j)
+            w = math.cast_like(w, p)
+            one = math.cast_like(one, p)
+            z = math.cast_like(z, p)
+
+        w2 = w**2
+
+        # The matrices are explicitly written, not generated to ensure PyTorch differentiation.
+        depolarizing_mats = [
+            [[one, z, z], [z, w, z], [z, z, w2]],
+            [[one, z, z], [z, w2, z], [z, z, w]],
+            [[z, one, z], [z, z, one], [one, z, z]],
+            [[z, w, z], [z, z, w2], [one, z, z]],
+            [[z, w2, z], [z, z, w], [one, z, z]],
+            [[z, z, one], [one, z, z], [z, one, z]],
+            [[z, z, w2], [one, z, z], [z, w, z]],
+            [[z, z, w], [one, z, z], [z, w2, z]],
+        ]
+
+        normalization = math.sqrt(p / 8 + math.eps)
+        Ks = [normalization * math.array(m, like=interface) for m in depolarizing_mats]
+        identity = math.sqrt(1 - p + math.eps) * math.array(
+            math.eye(QUDIT_DIM, dtype=complex), like=interface
+        )
+
+        return [identity] + Ks

tests/devices/qutrit_mixed/test_qutrit_mixed_preprocessing.py

@@ -125,9 +125,10 @@ def test_measurement_is_swapped_out(self, mp_fn, mp_cls, shots):
             (qml.GellMann(0, 1), False),
             (qml.Snapshot(), True),
             (qml.TRX(1.1, 0), True),
+            (qml.QutritDepolarizingChannel(0.4, 0), True),
         ],
     )
-    def test_accepted_observables(self, op, expected):
+    def test_accepted_operator(self, op, expected):
         """Test that stopping_condition works correctly"""
         res = stopping_condition(op)
         assert res == expected
@@ -136,6 +137,7 @@ def test_accepted_observables(self, op, expected):
         "obs, expected",
         [
             (qml.TShift(0), False),
+            (qml.QutritDepolarizingChannel(0.4, 0), False),
             (qml.GellMann(0, 1), True),
             (qml.Snapshot(), False),
             (qml.operation.Tensor(qml.GellMann(0, 1), qml.GellMann(3, 3)), True),
@@ -144,7 +146,7 @@ def test_accepted_observables(self, op, expected):
             (qml.ops.op_math.Prod(qml.GellMann(0, 1), qml.GellMann(3, 3)), True),
         ],
     )
-    def test_accepted_operator(self, obs, expected):
+    def test_accepted_observable(self, obs, expected):
         """Test that observable_stopping_condition works correctly"""
         res = observable_stopping_condition(obs)
         assert res == expected