From 332d09007f885f1e8ae06fbe99b7112d45f840d9 Mon Sep 17 00:00:00 2001 From: Gabriel Bottrill <78718539+Gabriel-Bottrill@users.noreply.github.com> Date: Wed, 15 May 2024 14:32:32 -0700 Subject: [PATCH] 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 Co-authored-by: Olivia Di Matteo <2068515+glassnotes@users.noreply.github.com> --- doc/introduction/operations.rst | 15 ++ doc/releases/changelog-dev.md | 6 + pennylane/devices/default_qutrit_mixed.py | 3 +- pennylane/ops/__init__.py | 1 + pennylane/ops/qutrit/__init__.py | 6 +- pennylane/ops/qutrit/channel.py | 224 ++++++++++++++++++ .../test_qutrit_mixed_preprocessing.py | 6 +- tests/ops/qutrit/test_qutrit_channel_ops.py | 174 ++++++++++++++ 8 files changed, 430 insertions(+), 5 deletions(-) create mode 100644 pennylane/ops/qutrit/channel.py create mode 100644 tests/ops/qutrit/test_qutrit_channel_ops.py diff --git a/doc/introduction/operations.rst b/doc/introduction/operations.rst index 9666fc58184..970e7b4ecd9 100644 --- a/doc/introduction/operations.rst +++ b/doc/introduction/operations.rst @@ -495,6 +495,21 @@ Qutrit State preparation :html:`` +.. _intro_ref_ops_qutrit_channels: + +Qutrit noisy channels +^^^^^^^^^^^^^^^^^^^^^^^^ + + +:html:`
` + +.. autosummary:: + :nosignatures: + + ~pennylane.QutritDepolarizingChannel + +:html:`
` + .. _intro_ref_ops_qutrit_obs: Qutrit Observables diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md index 0b6587eacb6..01db3008d03 100644 --- a/doc/releases/changelog-dev.md +++ b/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) +

Community contributions 🥳

+ +* ``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) +

Deprecations 👋

* ``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, diff --git a/pennylane/devices/default_qutrit_mixed.py b/pennylane/devices/default_qutrit_mixed.py index 7f98a9497c0..7c264a9d785 100644 --- a/pennylane/devices/default_qutrit_mixed.py +++ b/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", diff --git a/pennylane/ops/__init__.py b/pennylane/ops/__init__.py index ee1a9eefa3c..6beae6b802e 100644 --- a/pennylane/ops/__init__.py +++ b/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__) diff --git a/pennylane/ops/qutrit/__init__.py b/pennylane/ops/qutrit/__init__.py index 36aa7f22789..670f950ba45 100644 --- a/pennylane/ops/qutrit/__init__.py +++ b/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__) diff --git a/pennylane/ops/qutrit/channel.py b/pennylane/ops/qutrit/channel.py new file mode 100644 index 00000000000..8443fd2de44 --- /dev/null +++ b/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 `_] (Eq. 5). + The Kraus Matrices we use are adapted from [`2 `_] (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 diff --git a/tests/devices/qutrit_mixed/test_qutrit_mixed_preprocessing.py b/tests/devices/qutrit_mixed/test_qutrit_mixed_preprocessing.py index 27bf1331d17..9574fdbec25 100644 --- a/tests/devices/qutrit_mixed/test_qutrit_mixed_preprocessing.py +++ b/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 diff --git a/tests/ops/qutrit/test_qutrit_channel_ops.py b/tests/ops/qutrit/test_qutrit_channel_ops.py new file mode 100644 index 00000000000..be1025d0c26 --- /dev/null +++ b/tests/ops/qutrit/test_qutrit_channel_ops.py @@ -0,0 +1,174 @@ +# Copyright 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. +""" +Unit tests for the available built-in qutrit quantum channels. +""" +import numpy as np +import pytest +from numpy.linalg import matrix_power + +import pennylane as qml +from pennylane import numpy as pnp +from pennylane.ops.qutrit import channel + +QUDIT_DIM = 3 + + +class TestQutritDepolarizingChannel: + """Tests for the qutrit quantum channel QutritDepolarizingChannel""" + + @staticmethod + def get_expected_kraus_matrices(p): + """Gets the expected Kraus matrices given probability p.""" + w = np.exp(2j * np.pi / 3) + + X = np.array([[0, 1, 0], [0, 0, 1], [1, 0, 0]]) + Z = np.diag([1, w, w**2]) + + Ks = [np.sqrt(1 - p) * np.eye(QUDIT_DIM)] + for i in range(3): + for j in range(3): + if i == 0 and j == 0: + continue + Ks.append(np.sqrt(p / 8) * matrix_power(X, i) @ matrix_power(Z, j)) + return Ks + + def test_p_zero(self, tol): + """Test p=0 gives correct Kraus matrices""" + op = qml.QutritDepolarizingChannel + kraus_matrices = op(0, wires=0).kraus_matrices() + + assert len(kraus_matrices) == 9 + assert np.allclose(kraus_matrices[0], np.eye(QUDIT_DIM), atol=tol, rtol=0) + assert np.allclose(np.array(kraus_matrices[1:]), 0, atol=tol, rtol=0) + + def test_p_arbitrary(self, tol): + """Test p=0.1 gives correct Kraus matrices""" + p = 0.1 + kraus_matrices = qml.QutritDepolarizingChannel(p, wires=0).kraus_matrices() + expected_matrices = self.get_expected_kraus_matrices(p) + for kraus_matrix, expected_matrix in zip(kraus_matrices, expected_matrices): + assert np.allclose(kraus_matrix, expected_matrix, atol=tol, rtol=0) + + def test_p_invalid_parameter(self): + """Test that error is raised given an inappropriate p value.""" + with pytest.raises(ValueError, match="p must be in the interval"): + qml.QutritDepolarizingChannel(1.5, wires=0).kraus_matrices() + + @pytest.mark.parametrize("angle", np.linspace(0, 2 * np.pi, 7)) + def test_grad_depolarizing(self, angle): + """Test that analytical gradient is computed correctly for different states. Channel + grad recipes are independent of channel parameter""" + + dev = qml.device("default.qutrit.mixed") + prob = pnp.array(0.5, requires_grad=True) + + @qml.qnode(dev) + def circuit(p): + qml.TRX(angle, wires=0, subspace=(0, 1)) + qml.TRX(angle, wires=0, subspace=(1, 2)) + qml.QutritDepolarizingChannel(p, wires=0) + return qml.expval(qml.GellMann(0, 3) + qml.GellMann(0, 8)) + + expected_errorless = ( + (np.sqrt(3) - 3) * (1 - np.cos(2 * angle)) / 24 + - 2 / np.sqrt(3) * np.sin(angle / 2) ** 4 + + (np.sqrt(1 / 3) + 1) * np.cos(angle / 2) ** 2 + ) + + assert np.allclose(circuit(prob), ((prob - (1 / 9)) / (8 / 9)) * expected_errorless) + + gradient = np.squeeze(qml.grad(circuit)(prob)) + assert np.allclose(gradient, circuit(1) - circuit(0)) + assert np.allclose(gradient, -(9 / 8) * expected_errorless) + + @staticmethod + def expected_jac_fn(p): + """Gets the expected Jacobian of Kraus matrices given probability p.""" + w = np.exp(2j * np.pi / 3) + + X = np.array([[0, 1, 0], [0, 0, 1], [1, 0, 0]]) + Z = np.diag([1, w, w**2]) + + jacs = [-1 / (2 * np.sqrt(1 - p)) * np.eye(QUDIT_DIM)] + for i in range(3): + for j in range(3): + if i == 0 and j == 0: + continue + jacs.append(np.sqrt(1 / (32 * p)) * matrix_power(X, i) @ matrix_power(Z, j)) + + return jacs + + @staticmethod + def kraus_fn(p): + """Gets a matrix of the Kraus matrices to be tested.""" + return qml.math.stack(channel.QutritDepolarizingChannel(p, wires=0).kraus_matrices()) + + @staticmethod + def kraus_fn_real(p): + """Gets a matrix of the real part of the Kraus matrices to be tested.""" + return qml.math.real( + qml.math.stack(channel.QutritDepolarizingChannel(p, wires=0).kraus_matrices()) + ) + + @staticmethod + def kraus_fn_imag(p): + """Gets a matrix of the imaginary part of the Kraus matrices to be tested.""" + return qml.math.imag( + qml.math.stack(channel.QutritDepolarizingChannel(p, wires=0).kraus_matrices()) + ) + + @pytest.mark.autograd + def test_kraus_jac_autograd(self): + """Tests Jacobian of Kraus matrices using autograd.""" + p = pnp.array(0.43, requires_grad=True) + jac = qml.jacobian(self.kraus_fn_real)(p) + 1j * qml.jacobian(self.kraus_fn_imag)(p) + assert qml.math.allclose(jac, self.expected_jac_fn(p)) + + @pytest.mark.torch + def test_kraus_jac_torch(self): + """Tests Jacobian of Kraus matrices using torch.""" + import torch + + p = torch.tensor(0.43, requires_grad=True) + jacobian = torch.autograd.functional.jacobian + jac = jacobian(self.kraus_fn_real, p) + 1j * jacobian(self.kraus_fn_imag, p) + assert qml.math.allclose(jac, self.expected_jac_fn(p.detach().numpy())) + + @pytest.mark.tf + def test_kraus_jac_tf(self): + """Tests Jacobian of Kraus matrices using tensorflow.""" + import tensorflow as tf + + p = tf.Variable(0.43) + with tf.GradientTape() as real_tape: + real_out = self.kraus_fn_real(p) + with tf.GradientTape() as imag_tape: + imag_out = self.kraus_fn_imag(p) + + real_jac = qml.math.cast(real_tape.jacobian(real_out, p), complex) + imag_jac = qml.math.cast(imag_tape.jacobian(imag_out, p), complex) + jac = real_jac + 1j * imag_jac + assert qml.math.allclose(jac, self.expected_jac_fn(0.43)) + + @pytest.mark.jax + def test_kraus_jac_jax(self): + """Tests Jacobian of Kraus matrices using jax.""" + import jax + + jax.config.update("jax_enable_x64", True) + + p = jax.numpy.array(0.43, dtype=jax.numpy.complex128) + jac = jax.jacobian(self.kraus_fn, holomorphic=True)(p) + assert qml.math.allclose(jac, self.expected_jac_fn(p))