Qutrit channel depolarizing (#5502)

**Context:**
`default.qutrit.mixed` device has been added, but no channels have been
added so far, this adds the first channel to the device so that it may
be used for it's intended purpose of simulating noise.

**Description of the Change:**
Adds new channel module to `qml.ops.qutrit` package. Adding the first
qutrit channel `QutritDepolarizingChannel`.

**Benefits:**
Allows for `defualt.qutrit.mixed` to simulate depolarizing noise. Makes
`default.qutrit.mixed`, much more useful.

**Possible Drawbacks:**
N/A

**Related GitHub Issues:**
N/A

---------

Co-authored-by: Gabriel Bottrill <bottrill@student.ubc.ca>
Co-authored-by: Olivia Di Matteo <2068515+glassnotes@users.noreply.github.com>

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