Repository: PennyLaneAI/pennylane
Pull request #1766 (2021-10-19 00:50:15)

Add noisy channel: thermal relaxation
### Before submitting

Please complete the following checklist when submitting a PR:

- [x] All new features must include a unit test.
      If you've fixed a bug or added code that should be tested, add a test to the
      test directory!

- [x] All new functions and code must be clearly commented and documented.
      If you do make documentation changes, make sure that the docs build and
      render correctly by running `make docs`.

- [x] Ensure that the test suite passes, by running `make test`.

- [x] Add a new entry to the `doc/releases/changelog-dev.md` file, summarizing the
      change, and including a link back to the PR.

- [x] The PennyLane source code conforms to
      [PEP8 standards](https://www.python.org/dev/peps/pep-0008/).
      We check all of our code against [Pylint](https://www.pylint.org/).
      To lint modified files, simply `pip install pylint`, and then
      run `pylint pennylane/path/to/file.py`.

When all the above are checked, delete everything above the dashed
line and fill in the pull request template.

------------------------------------------------------------------------------------------------------------

**Context:**
Adding a new noisy channel which gives the thermal relaxation channel as part of the #971 issue. I got the information from the [supplementary](https://static-content.springer.com/esm/art%3A10.1038%2Fs41534-020-0272-6/MediaObjects/41534_2020_272_MOESM1_ESM.pdf) of [Quantum classifier with tailored quantum kernel](https://www.nature.com/articles/s41534-020-0272-6.pdf). I still need to figure out how to get the second case when T2 > T1.
**Description of the Change:**

**Benefits:**

**Possible Drawbacks:**

**Related GitHub Issues:**


josh146: Thanks @alejomonbar for this nice contribution! 🙂 Let us know if you have any questions, or when this PR is ready for review, and someone from the team will be on hand to help.
alejomonbar: > @alejomonbar I checked the different documents about thermal relaxation. As PennyLane has no Choi representation of quantum channel, we indeed need to construct the kraus matrices from the Choi matrix given in the paper. Therefore for the case where dephasing is more prevalent than relaxation (T2>T1), you will need indeed need this transformation. For the moment you can create this function directly in `channel.py` and work there. I will check if a better location is needed.
> 
> Thank you again for your contribution and feel free to ask any question. Let me know when the PR will be ready for review!

Thank you, I will look at that and come back when I implement it!
alejomonbar: > Thank you for the contribution, before I can give a full review you have different points to work on:
> 
> 1. Codefactor: in the folder `channel.py` minor things like trailinng white spaces and too many statements should be solved easily.
> 2. You can `pip install black` and run `black -l 100 your_file.py` on all the changed files.
> 3. For docstrings and documentation, there probably some small errors of formatting.
> 4. Some tests are failing, like `test_kraus_matrices_sum_identity`did you check it locally? Maybe there is small mistake in the Kraus operators or you need to adapt the test for the case. Same for the new tests are failing as well, do they pass locally for you?

Thank you @rmoyard , I tried it but I have a problem when trying to run the test I got the following message, `ModuleNotFoundError:No module named 'autoray'` and I installed it using `pip install autoray`, and looks like it is installed, but still, I got the same message.

rmoyard: > > Thank you for the contribution, before I can give a full review you have different points to work on:
> > 
> > 1. Codefactor: in the folder `channel.py` minor things like trailinng white spaces and too many statements should be solved easily.
> > 2. You can `pip install black` and run `black -l 100 your_file.py` on all the changed files.
> > 3. For docstrings and documentation, there probably some small errors of formatting.
> > 4. Some tests are failing, like `test_kraus_matrices_sum_identity`did you check it locally? Maybe there is small mistake in the Kraus operators or you need to adapt the test for the case. Same for the new tests are failing as well, do they pass locally for you?
> 
> Thank you @rmoyard , I tried it but I have a problem when trying to run the test I got the following message, `ModuleNotFoundError:No module named 'autoray'` and I installed it using `pip install autoray`, and looks like it is installed, but still, I got the same message.

Does it appear in the environment if you write `pip list`?
alejomonbar: > > > Thank you for the contribution, before I can give a full review you have different points to work on:
> > > 
> > > 1. Codefactor: in the folder `channel.py` minor things like trailinng white spaces and too many statements should be solved easily.
> > > 2. You can `pip install black` and run `black -l 100 your_file.py` on all the changed files.
> > > 3. For docstrings and documentation, there probably some small errors of formatting.
> > > 4. Some tests are failing, like `test_kraus_matrices_sum_identity`did you check it locally? Maybe there is small mistake in the Kraus operators or you need to adapt the test for the case. Same for the new tests are failing as well, do they pass locally for you?
> > 
> > 
> > Thank you @rmoyard , I tried it but I have a problem when trying to run the test I got the following message, `ModuleNotFoundError:No module named 'autoray'` and I installed it using `pip install autoray`, and looks like it is installed, but still, I got the same message.
> 
> Does it appear in the environment if you write `pip list`?

Yes, it does with version `autoray 0.2.5`
rmoyard: @alejomonbar Could you give me your system information? `import pennylane as qml; qml.about()`

Patch:
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index 75e0884e0..757ac23ff 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -3,6 +3,9 @@
 # Release 0.20.0-dev (development release)
 
 <h3>New features since last release</h3>
+* A thermal relaxation channel is added to the Noisy channels. The channel description can be 
+  found on the supplementary information of [Quantum classifier with tailored quantum kernels](https://arxiv.org/abs/1909.02611).
+  [(#1766)](https://github.com/PennyLaneAI/pennylane/pull/1766)
 
 <h3>Improvements</h3>
 
@@ -22,4 +25,4 @@
 
 This release contains contributions from (in alphabetical order):
 
-Christina Lee
+Christina Lee, Alejandro Montanez
diff --git a/pennylane/devices/default_mixed.py b/pennylane/devices/default_mixed.py
index 81bc9c9b5..db3e7d7e8 100644
--- a/pennylane/devices/default_mixed.py
+++ b/pennylane/devices/default_mixed.py
@@ -103,6 +103,7 @@ class DefaultMixed(QubitDevice):
         "QubitSum",
         "OrbitalRotation",
         "QFT",
+        "ThermalRelaxationError",
     }
 
     def __init__(self, wires, *, shots=None, cache=0, analytic=None):
diff --git a/pennylane/ops/channel.py b/pennylane/ops/channel.py
index 422568197..ab7ffd2c4 100644
--- a/pennylane/ops/channel.py
+++ b/pennylane/ops/channel.py
@@ -466,6 +466,147 @@ class QubitChannel(Channel):
         return K_list
 
 
+class ThermalRelaxationError(Channel):
+    r"""ThermalRelaxationError(pe, t1, t2, tg, wires)
+    Thermal relaxation error channel.
+
+    This channel is modelled by the following Kraus matrices:
+
+    Case :math:`T_2 \leq T_1`:
+
+    .. math::
+        K_0 = \sqrt{1 - p_z - p_{r0} - p_{r1}} \begin{bmatrix}
+                1 & 0 \\
+                0 & 1
+                \end{bmatrix}
+
+    .. math::
+        K_1 = \sqrt{p_z}\begin{bmatrix}
+                1 & 0  \\
+                0 & -1
+                \end{bmatrix}
+
+    .. math::
+        K_2 = \sqrt{p_{r0}}\begin{bmatrix}
+                1 & 0  \\
+                0 & 0
+                \end{bmatrix}
+
+    .. math::
+        K_3 = \sqrt{p_{r0}}\begin{bmatrix}
+                0 & 1  \\
+                0 & 0
+                \end{bmatrix}
+
+    .. math::
+        K_4 = \sqrt{p_{r1}}\begin{bmatrix}
+                0 & 0  \\
+                1 & 0
+                \end{bmatrix}
+
+    .. math::
+        K_5 = \sqrt{p_{r1}}\begin{bmatrix}
+                0 & 0  \\
+                0 & 1
+                \end{bmatrix}
+
+    where :math:`p_{r0} \in [0, 1]` is the probability of a reset to 0, :math:`p_{r1} \in [0, 1]` is the probability of
+    a reset to 1 error, :math:`p_z \in [0, 1]` is the probability of a phase flip (Pauli :math:`Z`) error.
+
+    Case :math:`T_2 > T_1`:
+    The Choi matrix is given by
+
+    .. math::
+        \Lambda = \begin{bmatrix}
+                        1 - p_e * p_{reset} & 0 & 0 & eT_2 \\
+                        0 & p_e * p_{reset} & 0 & 0 \\
+                        0 & 0 & (1 - p_e) * p_{reset} & 0 \\
+                        eT_2 & 0 & 0 & 1 - (1 - p_e) * p_{reset}
+                        \end{bmatrix}
+
+    .. math::
+        K_N = \sqrt{\lambda} \Phi(\nu_{\lambda})
+
+    where :math:`\lambda` are the eigenvalues of the Choi matrix, :math:`\nu_{\lambda}` are the eigenvectors of
+    the choi_matrix, and :math:`\Phi(x)` is a isomorphism from :math:`\mathbb{C}^{n^2}`
+    to :math:`\mathbb{C}^{n \times n}` with column-major order mapping.
+
+    **Details:**
+
+    * Number of wires: 1
+    * Number of parameters: 4
+
+    Args:
+        pe (float): exited state population.
+        t1 (float): the :math:`T_1` relaxation constant.
+        t2 (float): the :math:`T_2` dephasing constant.
+        tg (float): the gate time for relaxation error.
+        wires (Sequence[int] or int): the wire the channel acts on
+    """
+    num_params = 4
+    num_wires = 1
+    par_domain = "R"
+    grad_method = "F"
+
+    @classmethod
+    def _kraus_matrices(cls, *params):
+        pe = params[0]
+        t1 = params[1]
+        t2 = params[2]
+        tg = params[3]
+        if not 0.0 <= pe <= 1.0:
+            raise ValueError("pe must be between")
+        if tg < 0:
+            raise ValueError(f"Invalid gate_time ({tg} < 0)")
+        if t1 <= 0:
+            raise ValueError("Invalid T_1 relaxation time parameter: T_1 <= 0.")
+        if t2 <= 0:
+            raise ValueError("Invalid T_2 relaxation time parameter: T_2 <= 0.")
+        if t2 - 2 * t1 > 0:
+            raise ValueError("Invalid T_2 relaxation time parameter: T_2 greater than 2 * T_1.")
+        # T1 relaxation rate
+        eT1 = np.exp(-tg / t1)
+        p_reset = 1 - eT1
+        # T2 dephasing rate
+        eT2 = np.exp(-tg / t2)
+
+        if t2 <= t1:
+            pz = (1 - p_reset) * (1 - eT2 / eT1) / 2
+            pr0 = (1 - pe) * p_reset
+            pr1 = pe * p_reset
+            pid = 1 - pz - pr0 - pr1
+
+            K0 = np.sqrt(pid) * np.eye(2)
+            K1 = np.sqrt(pz) * np.array([[1, 0], [0, -1]])
+            K2 = np.sqrt(pr0) * np.array([[1, 0], [0, 0]])
+            K3 = np.sqrt(pr0) * np.array([[0, 1], [0, 0]])
+            K4 = np.sqrt(pr1) * np.array([[0, 0], [1, 0]])
+            K5 = np.sqrt(pr1) * np.array([[0, 0], [0, 1]])
+
+            K = [K0, K1, K2, K3, K4, K5]
+        else:
+            e0 = p_reset * pe
+            v0 = np.array([[0, 0], [1, 0]])
+            K0 = np.sqrt(e0) * v0
+            e1 = -p_reset * pe + p_reset
+            v1 = np.array([[0, 1], [0, 0]])
+            K1 = np.sqrt(e1) * v1
+            common_term = np.sqrt(
+                4 * eT2 ** 2 + 4 * p_reset ** 2 * pe ** 2 - 4 * p_reset ** 2 * pe + p_reset ** 2
+            )
+            e2 = 1 - p_reset / 2 - common_term / 2
+            term2 = 2 * eT2 / (2 * p_reset * pe - p_reset - common_term)
+            v2 = np.array([[term2, 0], [0, 1]]) / np.sqrt(term2 ** 2 + 1)
+            K2 = np.sqrt(e2) * v2
+            term3 = 2 * eT2 / (2 * p_reset * pe - p_reset + common_term)
+            e3 = 1 - p_reset / 2 + common_term / 2
+            v3 = np.array([[term3, 0], [0, 1]]) / np.sqrt(term3 ** 2 + 1)
+            K3 = np.sqrt(e3) * v3
+
+            K = [K0, K1, K2, K3]
+        return K
+
+
 __qubit_channels__ = {
     "AmplitudeDamping",
     "GeneralizedAmplitudeDamping",
@@ -475,6 +616,7 @@ __qubit_channels__ = {
     "PhaseFlip",
     "ResetError",
     "QubitChannel",
+    "ThermalRelaxationError",
 }
 
 __all__ = list(__qubit_channels__)


Test patch:
diff --git a/tests/ops/test_channel_ops.py b/tests/ops/test_channel_ops.py
index 022bb2d7f..f612287ce 100644
--- a/tests/ops/test_channel_ops.py
+++ b/tests/ops/test_channel_ops.py
@@ -14,8 +14,8 @@
 """
 Unit tests for the available built-in quantum channels.
 """
-import pytest
 import functools
+import pytest
 import numpy as np
 import pennylane as qml
 from pennylane.ops import channel
@@ -33,6 +33,7 @@ ch_list = [
     channel.PhaseFlip,
     channel.DepolarizingChannel,
     channel.ResetError,
+    channel.ThermalRelaxationError,
 ]
 
 
@@ -41,12 +42,15 @@ class TestChannels:
 
     @pytest.mark.parametrize("ops", ch_list)
     @pytest.mark.parametrize("p", [0, 0.1, 1])
-    def test_kraus_matrices_sum_identity(self, ops, p, tol):
+    @pytest.mark.parametrize("tr_args", [[100e-6, 100e-6, 20e-9], [100e-6, 120e-6, 20e-9]])
+    def test_kraus_matrices_sum_identity(self, ops, p, tr_args, tol):
         """Test channels are trace-preserving"""
         if ops.__name__ == "GeneralizedAmplitudeDamping":
             op = ops(p, p, wires=0)
         elif ops.__name__ == "ResetError":
             op = ops(p / 2, p / 3, wires=0)
+        elif ops.__name__ == "ThermalRelaxationError":
+            op = ops(p, *tr_args, wires=0)
         else:
             op = ops(p, wires=0)
         K_list = op.kraus_matrices
@@ -364,3 +368,166 @@ class TestQubitChannel:
         K_list2 = [np.sqrt(p) * Y, np.sqrt(1 - p) * np.eye(2)]
         with pytest.raises(ValueError, match="Only trace preserving channels can be applied."):
             channel.QubitChannel(K_list2 * 2, wires=0)
+
+
+class TestThermalRelaxationError:
+    """Tests for the quantum channel ThermalRelaxationError"""
+
+    @pytest.mark.parametrize(
+        "pe,t1,t2,tg",
+        list(
+            zip(
+                [0.2, 0.4, 0.6, 0.0],
+                [100e-6, 50e-6, 80e-6, np.inf],
+                [80e-6, 40e-6, 80e-6, 50e-6],
+                [20e-9, 40e-9, 40e-6, 40e-9],
+            )
+        ),
+    )
+    def test_t2_le_t1_arbitrary(self, pe, t1, t2, tg, tol):
+        """Test that various values of pe, t1, t2, and tg  for t2 <= t1 give correct Kraus matrices"""
+
+        op = channel.ThermalRelaxationError
+
+        eT1 = np.exp(-tg / t1)
+        p_reset = 1 - eT1
+        eT2 = np.exp(-tg / t2)
+        pz = (1 - p_reset) * (1 - eT2 / eT1) / 2
+        pr0 = (1 - pe) * p_reset
+        pr1 = pe * p_reset
+        pid = 1 - pz - pr0 - pr1
+
+        expected_K0 = np.sqrt(pid) * np.eye(2)
+        assert np.allclose(
+            op(pe, t1, t2, tg, wires=0).kraus_matrices[0], expected_K0, atol=tol, rtol=0
+        )
+
+        expected_K1 = np.sqrt(pz) * np.array([[1, 0], [0, -1]])
+        assert np.allclose(
+            op(pe, t1, t2, tg, wires=0).kraus_matrices[1], expected_K1, atol=tol, rtol=0
+        )
+
+        expected_K2 = np.sqrt(pr0) * np.array([[1, 0], [0, 0]])
+        assert np.allclose(
+            op(pe, t1, t2, tg, wires=0).kraus_matrices[2], expected_K2, atol=tol, rtol=0
+        )
+
+        expected_K3 = np.sqrt(pr0) * np.array([[0, 1], [0, 0]])
+        assert np.allclose(
+            op(pe, t1, t2, tg, wires=0).kraus_matrices[3], expected_K3, atol=tol, rtol=0
+        )
+
+        expected_K4 = np.sqrt(pr1) * np.array([[0, 0], [1, 0]])
+        assert np.allclose(
+            op(pe, t1, t2, tg, wires=0).kraus_matrices[4], expected_K4, atol=tol, rtol=0
+        )
+
+        expected_K5 = np.sqrt(pr1) * np.array([[0, 0], [0, 1]])
+        assert np.allclose(
+            op(pe, t1, t2, tg, wires=0).kraus_matrices[5], expected_K5, atol=tol, rtol=0
+        )
+
+    @pytest.mark.parametrize(
+        "pe,t1,t2,tg",
+        list(
+            zip(
+                [0.8, 0.5, 0.0],
+                [100e-6, 50e-6, 80e-6],
+                [120e-6, 100e-6, 90e-6],
+                [20e-9, 40e-9, 90e-6],
+            )
+        ),
+    )
+    def test_t2_g_t1_arbitrary(self, pe, t1, t2, tg, tol):
+        """Test that various values of pe, t1, t2, and tg  for t2 > t1 give correct Kraus matrices"""
+
+        op = channel.ThermalRelaxationError
+
+        if t1 == np.inf:
+            eT1 = 0
+            p_reset = 0
+        else:
+            eT1 = np.exp(-tg / t1)
+            p_reset = 1 - eT1
+        if t2 == np.inf:
+            eT2 = 1
+        else:
+            eT2 = np.exp(-tg / t2)
+
+        e0 = p_reset * pe
+        v0 = np.array([[0], [1], [0], [0]])
+        e1 = -p_reset * pe + p_reset
+        v1 = np.array([[0], [0], [1], [0]])
+        common_term = np.sqrt(
+            4 * eT2 ** 2 + 4 * p_reset ** 2 * pe ** 2 - 4 * p_reset ** 2 * pe + p_reset ** 2
+        )
+        e2 = 1 - p_reset / 2 - common_term / 2
+        term2 = 2 * eT2 / (2 * p_reset * pe - p_reset - common_term)
+        v2 = np.array([[term2], [0], [0], [1]]) / np.sqrt(term2 ** 2 + 1 ** 2)
+        term3 = 2 * eT2 / (2 * p_reset * pe - p_reset + common_term)
+        e3 = 1 - p_reset / 2 + common_term / 2
+        v3 = np.array([[term3], [0], [0], [1]]) / np.sqrt(term3 ** 2 + 1 ** 2)
+
+        expected_K0 = np.sqrt(e0) * v0.reshape(2, 2, order="F")
+        assert np.allclose(
+            op(pe, t1, t2, tg, wires=0).kraus_matrices[0], expected_K0, atol=tol, rtol=0
+        )
+
+        expected_K1 = np.sqrt(e1) * v1.reshape(2, 2, order="F")
+        assert np.allclose(
+            op(pe, t1, t2, tg, wires=0).kraus_matrices[1], expected_K1, atol=tol, rtol=0
+        )
+
+        expected_K2 = np.sqrt(e2) * v2.reshape(2, 2, order="F")
+        assert np.allclose(
+            op(pe, t1, t2, tg, wires=0).kraus_matrices[2], expected_K2, atol=tol, rtol=0
+        )
+
+        expected_K3 = np.sqrt(e3) * v3.reshape(2, 2, order="F")
+        assert np.allclose(
+            op(pe, t1, t2, tg, wires=0).kraus_matrices[3], expected_K3, atol=tol, rtol=0
+        )
+
+    def test_pe_invalid_parameter(self):
+        with pytest.raises(ValueError, match="pe must be between"):
+            channel.ThermalRelaxationError(1.5, 100e-6, 100e-6, 20e-9, wires=0).kraus_matrices
+
+    def test_T2_g_2T1_invalid_parameter(self):
+        with pytest.raises(ValueError, match="Invalid T_2 relaxation time parameter"):
+            channel.ThermalRelaxationError(0.3, 100e-6, np.inf, 20e-9, wires=0).kraus_matrices
+
+    def test_T1_le_0_invalid_parameter(self):
+        with pytest.raises(ValueError, match="Invalid T_1 relaxation time parameter"):
+            channel.ThermalRelaxationError(0.3, -50e-6, np.inf, 20e-9, wires=0).kraus_matrices
+
+    def test_T2_le_0_invalid_parameter(self):
+        with pytest.raises(ValueError, match="Invalid T_2 relaxation time parameter"):
+            channel.ThermalRelaxationError(0.3, 100e-6, 0, 20e-9, wires=0).kraus_matrices
+
+    def test_tg_le_0_invalid_parameter(self):
+        with pytest.raises(ValueError, match="Invalid gate_time"):
+            channel.ThermalRelaxationError(0.3, 100e-6, 100e-6, -20e-9, wires=0).kraus_matrices
+
+    @pytest.mark.parametrize("angle", np.linspace(0, 2 * np.pi, 7))
+    def test_grad_thermal_relaxation_error(self, angle, tol):
+        """Test that gradient is computed correctly for different states. Channel
+        grad recipes are independent of channel parameter"""
+
+        dev = qml.device("default.mixed", wires=1)
+        pe = 0.0
+
+        @qml.qnode(dev)
+        def circuit(pe):
+            qml.RX(angle, wires=0)
+            qml.ThermalRelaxationError(pe, 120e-6, 100e-6, 20e-9, wires=0)
+            return qml.expval(qml.PauliZ(0))
+
+        gradient = np.squeeze(qml.grad(circuit)(pe))
+        assert np.allclose(
+            gradient,
+            np.array(
+                [
+                    (1 / 0.1) * (circuit(0.1) - circuit(0.0)),
+                ]
+            ),
+        )