Gaussian random noise channel (#511)

**Context:**
Added the Gaussian random noise channel.

**Description of the Change:**
The main object is added as a `Channel`.

**Benefits:**
Clear.

**Possible Drawbacks:**
None.

**Related GitHub Issues:**
None

mrmustard/lab_dev/transformations/__init__.py

@@ -20,11 +20,12 @@
 from .attenuator import *
 from .base import *
 from .bsgate import *
+from .cft import *
 from .dgate import *
 from .fockdamping import *
 from .ggate import *
+from .gaussrandnoise import *
 from .identity import *
 from .rgate import *
 from .s2gate import *
 from .sgate import *
-from .cft import *

mrmustard/lab_dev/transformations/gaussrandnoise.py

@@ -0,0 +1,82 @@
+# 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.
+
+"""
+The class representing a Gaussian random noise channel.
+"""
+
+from __future__ import annotations
+from typing import Sequence
+from mrmustard import math, settings
+from mrmustard.utils.typing import RealMatrix
+from .base import Channel
+from ...physics.representations import Bargmann
+from ...physics import triples
+from ..utils import make_parameter
+
+__all__ = ["GaussRandNoise"]
+
+
+class GaussRandNoise(Channel):
+    r"""
+    The Gaussian random noise channel.
+
+    The number of modes must match half of the size of the Y matrix.
+    .. code-block ::
+
+        >>> import numpy as np
+        >>> from mrmustard.lab_dev import GaussRandNoise
+
+        >>> channel = GaussRandNoise(modes=[1, 2], Y = .2 * np.eye(4))
+        >>> assert channel.modes == [1, 2]
+        >>> assert np.allclose(channel.Y.value, .2 * np.eye(4))
+
+    Args:
+        modes: The modes the channel is applied to
+        Y: The Y matrix of the Gaussian random noise
+        Y_train: whether the Y matrix is a trainable variable
+
+    ..details::
+        The Bargmann representation of the channel is computed via the formulas provided in the paper:
+        https://arxiv.org/pdf/2209.06069
+
+        The channel maps an inout covariance matrix ``cov`` as
+        ..math::
+                cov \mapsto cov + Y.
+    """
+
+    short_name = "GRN"
+
+    def __init__(
+        self,
+        modes: Sequence[int],
+        Y: RealMatrix,
+        Y_trainable: bool = False,
+    ):
+
+        if Y.shape[-1] // 2 != len(modes):
+            raise ValueError(
+                f"The number of modes {len(modes)} does not match the dimension of the "
+                f"Y matrix {Y.shape[-1] // 2}."
+            )
+
+        if (math.real(math.eigvals(Y)) >= -settings.ATOL).min() == 0:
+            raise ValueError("The input Y matrix has negative eigen-values.")
+
+        super().__init__(modes_out=modes, modes_in=modes, name="GRN")
+        self._add_parameter(make_parameter(Y_trainable, value=Y, name="Y", bounds=(None, None)))
+
+        self._representation = Bargmann.from_function(
+            fn=triples.gaussian_random_noise_Abc, Y=self.Y.value
+        )

mrmustard/physics/triples.py

@@ -22,7 +22,7 @@
 import numpy as np
 
 from mrmustard import math, settings
-from mrmustard.utils.typing import Matrix, Vector, Scalar
+from mrmustard.utils.typing import Matrix, Vector, Scalar, RealMatrix
 from mrmustard.physics.gaussian_integrals import complex_gaussian_integral_2
 
 
@@ -605,6 +605,58 @@ def fock_damping_Abc(
     return A, b, c
 
 
+def gaussian_random_noise_Abc(Y: RealMatrix) -> Union[Matrix, Vector, Scalar]:
+    r"""
+    The triple (A, b, c) for the gaussian random noise channel.
+
+    Args:
+        Y: the Y matrix of the Gaussian random noise channel.
+
+    Returns:
+        The ``(A, b, c)`` triple of the Gaussian random noise channel.
+    """
+    m = Y.shape[-1] // 2
+    xi = math.eye(2 * m, dtype=math.complex128) + Y / settings.HBAR
+    xi_inv = math.inv(xi)
+    xi_inv_in_blocks = math.block(
+        [[math.eye(2 * m) - xi_inv, xi_inv], [xi_inv, math.eye(2 * m) - xi_inv]]
+    )
+    R = (
+        1
+        / math.sqrt(complex(2))
+        * math.block(
+            [
+                [
+                    math.eye(m, dtype=math.complex128),
+                    1j * math.eye(m, dtype=math.complex128),
+                    math.zeros((m, 2 * m), dtype=math.complex128),
+                ],
+                [
+                    math.zeros((m, 2 * m), dtype=math.complex128),
+                    math.eye(m, dtype=math.complex128),
+                    -1j * math.eye(m, dtype=math.complex128),
+                ],
+                [
+                    math.eye(m, dtype=math.complex128),
+                    -1j * math.eye(m, dtype=math.complex128),
+                    math.zeros((m, 2 * m), dtype=math.complex128),
+                ],
+                [
+                    math.zeros((m, 2 * m), dtype=math.complex128),
+                    math.eye(m, dtype=math.complex128),
+                    1j * math.eye(m, dtype=math.complex128),
+                ],
+            ]
+        )
+    )
+
+    A = math.Xmat(2 * m) @ R @ xi_inv_in_blocks @ math.conj(R).T
+    b = math.zeros(4 * m)
+    c = 1 / math.sqrt(math.det(xi))
+
+    return A, b, c
+
+
 def bargmann_to_quadrature_Abc(n_modes: int, phi: float) -> tuple[Matrix, Vector, Scalar]:
     r"""
     The ``(A, b, c)`` triple of the multi-mode kernel :math:`\langle \vec{p}|\vec{z} \rangle` between bargmann representation with ABC Ansatz form and quadrature representation with ABC Ansatz.

tests/test_lab_dev/test_transformations/test_gaussrandnoise.py

@@ -0,0 +1,49 @@
+# 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.
+
+"""Tests for the ``GaussRandNoise`` class."""
+
+# pylint: disable=protected-access, missing-function-docstring, expression-not-assigned
+
+import numpy as np
+
+from mrmustard import math
+from mrmustard.lab_dev.states import DM
+from mrmustard.lab_dev.transformations import GaussRandNoise
+
+
+class TestGRN:
+    r"""
+    Tests for the ``GaussRandNoise`` class.
+    """
+
+    def test_init(self):
+        "Tests the GaussRandNoise initialization."
+
+        a = np.random.random((2, 2))
+        grn = GaussRandNoise([0], a @ a.T)
+        assert grn.name == "GRN"
+        assert grn.modes == [0]
+
+    def test_grn(self):
+        "Tests if the A matrix of GaussRandNoise is computed correctly."
+        a = np.random.random((4, 4))
+        Y = a @ a.T
+        phi = GaussRandNoise([0, 1], Y)
+
+        _, Y_ans = phi.XY
+
+        assert math.allclose(Y_ans, Y)
+        assert phi.is_physical
+        assert math.allclose((DM.random([0, 1]) >> phi).probability, 1.0)

tests/test_physics/test_triples.py

@@ -336,3 +336,21 @@ def test_attenuator_kraus_Abc(self, eta):
         B = Bargmann(*triples.attenuator_kraus_Abc(eta))
         Att = Bargmann(*triples.attenuator_Abc(eta))
         assert B[2] @ B[2] == Att
+
+    def test_gaussian_random_noise_Abc(self):
+
+        A, b, c = triples.gaussian_random_noise_Abc(np.eye(2))
+        A_by_hand = np.array(
+            [
+                [0.0, 0.5, 0.5, 0.0],
+                [0.5, 0.0, 0.0, 0.5],
+                [0.5, 0.0, 0.0, 0.5],
+                [0.0, 0.5, 0.5, 0.0],
+            ]
+        )
+        b_by_hand = np.zeros(4)
+        c_by_hand = 0.5
+
+        assert math.allclose(A, A_by_hand)
+        assert math.allclose(b, b_by_hand)
+        assert math.allclose(c, c_by_hand)