diff --git a/.github/CHANGELOG.md b/.github/CHANGELOG.md
index 66a5f92fb..b0baf5836 100644
--- a/.github/CHANGELOG.md
+++ b/.github/CHANGELOG.md
@@ -5,20 +5,6 @@
   [(#128)](https://github.com/XanaduAI/MrMustard/issues/128)
 
 * States in Gaussian and Fock representation now can be concatenated.
-  [(#130)](https://github.com/XanaduAI/MrMustard/pull/130)
-
-* Parameter passthrough allows to use custom parameters in the model, that is, objects accept correlated parameters. For example, 
-    ```python
-    from mrmustard.lab.gates import Sgate, BSgate
-    
-    BS = BSgate(theta=np.pi/4, theta_trainable=True)[0,1]
-    S0 = Sgate(r=BS.theta)[0]
-    S1 = Sgate(r=-BS.theta)[1]
-    
-    circ = S0 >> S1 >> BS
-    ```
-  [(#131)](https://github.com/XanaduAI/MrMustard/pull/131)
-
   ```python
   from mrmustard.lab.states import Gaussian, Fock'
   from mrmustard.lab.gates import Attenuator
@@ -35,6 +21,22 @@
 
   mixed_state.dm()
   ```
+  [(#130)](https://github.com/XanaduAI/MrMustard/pull/130)
+
+* Parameter passthrough allows to use custom parameters in the model, that is, objects accept correlated parameters. For example, 
+    ```python
+    from mrmustard.lab.gates import Sgate, BSgate
+    
+    BS = BSgate(theta=np.pi/4, theta_trainable=True)[0,1]
+    S0 = Sgate(r=BS.theta)[0]
+    S1 = Sgate(r=-BS.theta)[1]
+    
+    circ = S0 >> S1 >> BS
+    ```
+  [(#131)](https://github.com/XanaduAI/MrMustard/pull/131)
+
+* Adds the new trainable gate `RealInterferometer`: an interferometer that doesn't mix the q and p quadratures
+  [(#132)](https://github.com/XanaduAI/MrMustard/pull/132)
 
 ### Breaking changes
 
diff --git a/mrmustard/lab/gates.py b/mrmustard/lab/gates.py
index 5bea0c2fe..f6cb2bc7f 100644
--- a/mrmustard/lab/gates.py
+++ b/mrmustard/lab/gates.py
@@ -23,9 +23,12 @@
 from mrmustard import settings
 from mrmustard.lab.abstract import Transformation
 from mrmustard.utils.parametrized import Parametrized
-from mrmustard.utils import training
 from mrmustard.physics import gaussian
 
+from mrmustard.math import Math
+
+math = Math()
+
 __all__ = [
     "Dgate",
     "Sgate",
@@ -38,6 +41,7 @@
     "CZgate",
     "CXgate",
     "Interferometer",
+    "RealInterferometer",
     "Attenuator",
     "Amplifier",
     "AdditiveNoise",
@@ -429,7 +433,8 @@ class Interferometer(Parametrized, Transformation):
     It corresponds to a Ggate with zero mean and a ``2N x 2N`` orthogonal symplectic matrix.
 
     Args:
-        orthogonal (2d array): a valid orthogonal matrix. For N modes it must have shape `(2N,2N)`
+        orthogonal (2d array, optional): a valid orthogonal matrix. For N modes it must have shape `(2N,2N)`.
+            If set to `None` a random orthogonal matrix is used.
         orthogonal_trainable (bool): whether orthogonal is a trainable variable
     """
 
@@ -440,7 +445,8 @@ def __init__(
         orthogonal_trainable: bool = False,
     ):
         if orthogonal is None:
-            orthogonal = training.new_orthogonal(num_modes=num_modes)
+            U = math.random_unitary(num_modes)
+            orthogonal = math.block([[math.real(U), -math.imag(U)], [math.imag(U), math.real(U)]])
         super().__init__(
             orthogonal=orthogonal,
             orthogonal_trainable=orthogonal_trainable,
@@ -454,9 +460,59 @@ def X_matrix(self):
         return self.orthogonal
 
     def _validate_modes(self, modes):
-        if len(modes) != self.orthogonal.shape[1] // 2:
+        if len(modes) != self.orthogonal.shape[-1] // 2:
+            raise ValueError(
+                f"Invalid number of modes: {len(modes)} (should be {self.orthogonal.shape[-1] // 2})"
+            )
+
+    @property
+    def trainable_parameters(self) -> Dict[str, List[Trainable]]:
+        return {
+            "symplectic": [],
+            "orthogonal": [self.orthogonal] if self._orthogonal_trainable else [],
+            "euclidean": [],
+        }
+
+
+class RealInterferometer(Parametrized, Transformation):
+    r"""N-mode interferometer with a real unitary matrix (or block-diagonal orthogonal matrix).
+    Does not mix q's and p's.
+
+    Args:
+        orthogonal (2d array, optional): a valid orthogonal matrix. For N modes it must have shape `(N,N)`.
+            If set to `None` a random orthogonal matrix is used.
+        orthogonal_trainable (bool): whether orthogonal is a trainable variable
+    """
+
+    def __init__(
+        self,
+        num_modes: int,
+        orthogonal: Optional[Tensor] = None,
+        orthogonal_trainable: bool = False,
+    ):
+        if orthogonal is None:
+            orthogonal = math.random_orthogonal(num_modes)
+        super().__init__(
+            orthogonal=orthogonal,
+            orthogonal_trainable=orthogonal_trainable,
+            orthogonal_bounds=(None, None),
+            modes=list(range(num_modes)),
+        )
+        self.is_gaussian = True
+
+    @property
+    def X_matrix(self):
+        return math.block(
+            [
+                [self.orthogonal, math.zeros_like(self.orthogonal)],
+                [math.zeros_like(self.orthogonal), self.orthogonal],
+            ]
+        )
+
+    def _validate_modes(self, modes):
+        if len(modes) != self.orthogonal.shape[-1]:
             raise ValueError(
-                f"Invalid number of modes: {len(modes)} (should be {self.orthogonal.shape[1] // 2})"
+                f"Invalid number of modes: {len(modes)} (should be {self.orthogonal.shape[-1]})"
             )
 
     @property
@@ -487,7 +543,7 @@ def __init__(
         symplectic_trainable: bool = False,
     ):
         if symplectic is None:
-            symplectic = training.new_symplectic(num_modes=num_modes)
+            symplectic = math.random_symplectic(num_modes)
         super().__init__(
             symplectic=symplectic,
             symplectic_trainable=symplectic_trainable,
diff --git a/mrmustard/lab/states.py b/mrmustard/lab/states.py
index 1ed27af22..9635552df 100644
--- a/mrmustard/lab/states.py
+++ b/mrmustard/lab/states.py
@@ -22,7 +22,6 @@
 from mrmustard.lab.abstract import State
 from mrmustard.physics import gaussian, fock
 from mrmustard.utils.parametrized import Parametrized
-from mrmustard.utils import training
 from mrmustard.math import Math
 
 math = Math()
@@ -421,7 +420,7 @@ def __init__(
         normalize: bool = False,
     ):
         if symplectic is None:
-            symplectic = training.new_symplectic(num_modes=num_modes)
+            symplectic = math.random_symplectic(num_modes=num_modes)
         if eigenvalues is None:
             eigenvalues = gaussian.math.ones(num_modes) * settings.HBAR / 2
         if math.any(math.atleast_1d(eigenvalues) < settings.HBAR / 2):
diff --git a/mrmustard/math/math_interface.py b/mrmustard/math/math_interface.py
index 2f45d3732..9e994e9a7 100644
--- a/mrmustard/math/math_interface.py
+++ b/mrmustard/math/math_interface.py
@@ -19,7 +19,7 @@
 from itertools import product
 import numpy as np
 from scipy.special import binom
-from scipy.stats import unitary_group
+from scipy.stats import unitary_group, ortho_group
 
 from mrmustard.types import (
     List,
@@ -68,7 +68,6 @@ def abs(self, array: Tensor) -> Tensor:
         Returns:
             array: absolute value of array
         """
-        ...
 
     @abstractmethod
     def any(self, array: Tensor) -> bool:
@@ -80,7 +79,6 @@ def any(self, array: Tensor) -> bool:
         Returns:
             bool: True if any element of array is True
         """
-        ...
 
     @abstractmethod
     def arange(self, start: int, limit: int = None, delta: int = 1, dtype: Any = None) -> Tensor:
@@ -95,7 +93,7 @@ def arange(self, start: int, limit: int = None, delta: int = 1, dtype: Any = Non
         Returns:
             array: array of evenly spaced values
         """
-        ...  # NOTE: is float64 by default
+        # NOTE: is float64 by default
 
     @abstractmethod
     def asnumpy(self, tensor: Tensor) -> Tensor:
@@ -107,7 +105,6 @@ def asnumpy(self, tensor: Tensor) -> Tensor:
         Returns:
             array: NumPy array
         """
-        ...
 
     @abstractmethod
     def assign(self, tensor: Tensor, value: Tensor) -> Tensor:
@@ -120,7 +117,6 @@ def assign(self, tensor: Tensor, value: Tensor) -> Tensor:
         Returns:
             array: tensor with value assigned
         """
-        ...
 
     @abstractmethod
     def astensor(self, array: Tensor, dtype: str) -> Tensor:
@@ -133,7 +129,6 @@ def astensor(self, array: Tensor, dtype: str) -> Tensor:
         Returns:
             array: tensor with dtype
         """
-        ...
 
     @abstractmethod
     def atleast_1d(self, array: Tensor, dtype=None) -> Tensor:
@@ -146,7 +141,6 @@ def atleast_1d(self, array: Tensor, dtype=None) -> Tensor:
         Returns:
             array: array with at least one dimension
         """
-        ...
 
     @abstractmethod
     def cast(self, array: Tensor, dtype) -> Tensor:
@@ -159,7 +153,6 @@ def cast(self, array: Tensor, dtype) -> Tensor:
         Returns:
             array: array cast to dtype
         """
-        ...
 
     @abstractmethod
     def clip(self, array: Tensor, a_min: float, a_max: float) -> Tensor:
@@ -173,7 +166,6 @@ def clip(self, array: Tensor, a_min: float, a_max: float) -> Tensor:
         Returns:
             array: clipped array
         """
-        ...
 
     @abstractmethod
     def concat(self, values: Sequence[Tensor], axis: int) -> Tensor:
@@ -186,7 +178,6 @@ def concat(self, values: Sequence[Tensor], axis: int) -> Tensor:
         Returns:
             array: concatenated values
         """
-        ...
 
     @abstractmethod
     def conj(self, array: Tensor) -> Tensor:
@@ -198,7 +189,6 @@ def conj(self, array: Tensor) -> Tensor:
         Returns:
             array: complex conjugate of array
         """
-        ...
 
     @abstractmethod
     def constraint_func(
@@ -219,7 +209,6 @@ def constraint_func(
         Returns:
             function: constraint function
         """
-        ...
 
     @abstractmethod
     def convolution(
@@ -240,7 +229,6 @@ def convolution(
         Returns:
             array: convolved array
         """
-        ...
 
     @abstractmethod
     def cos(self, array: Tensor) -> Tensor:
@@ -252,7 +240,6 @@ def cos(self, array: Tensor) -> Tensor:
         Returns:
             array: cosine of array
         """
-        ...
 
     @abstractmethod
     def cosh(self, array: Tensor) -> Tensor:
@@ -264,7 +251,6 @@ def cosh(self, array: Tensor) -> Tensor:
         Returns:
             array: hyperbolic cosine of array
         """
-        ...
 
     @abstractmethod
     def det(self, matrix: Tensor) -> Tensor:
@@ -276,7 +262,6 @@ def det(self, matrix: Tensor) -> Tensor:
         Returns:
             determinant of matrix
         """
-        ...
 
     @abstractmethod
     def diag(self, array: Tensor, k: int) -> Tensor:
@@ -289,7 +274,6 @@ def diag(self, array: Tensor, k: int) -> Tensor:
         Returns:
             array: array with array inserted into the kth diagonal
         """
-        ...
 
     @abstractmethod
     def diag_part(self, array: Tensor) -> Tensor:
@@ -301,7 +285,10 @@ def diag_part(self, array: Tensor) -> Tensor:
         Returns:
             array: array of the main diagonal of array
         """
-        ...
+
+    @abstractmethod
+    def eigvals(self, tensor: Tensor) -> Tensor:
+        r"""Returns the eigenvalues of a matrix."""
 
     @abstractmethod
     def einsum(self, string: str, *tensors) -> Tensor:
@@ -314,7 +301,6 @@ def einsum(self, string: str, *tensors) -> Tensor:
         Returns:
             array: result of the Einstein summation convention
         """
-        ...
 
     @abstractmethod
     def exp(self, array: Tensor) -> Tensor:
@@ -326,7 +312,6 @@ def exp(self, array: Tensor) -> Tensor:
         Returns:
             array: exponential of array
         """
-        ...
 
     @abstractmethod
     def expand_dims(self, array: Tensor, axis: int) -> Tensor:
@@ -339,7 +324,6 @@ def expand_dims(self, array: Tensor, axis: int) -> Tensor:
         Returns:
             array: array with an additional dimension inserted at the given axis
         """
-        ...
 
     @abstractmethod
     def expm(self, matrix: Tensor) -> Tensor:
@@ -351,7 +335,6 @@ def expm(self, matrix: Tensor) -> Tensor:
         Returns:
             matrix: exponential of matrix
         """
-        ...
 
     @abstractmethod
     def eye(self, size: int, dtype) -> Tensor:
@@ -364,12 +347,10 @@ def eye(self, size: int, dtype) -> Tensor:
         Returns:
             matrix: identity matrix
         """
-        ...
 
     @abstractmethod
     def from_backend(self, value: Any) -> bool:
         r"""Returns whether the given tensor is a tensor of the concrete backend."""
-        ...
 
     @abstractmethod
     def gather(self, array: Tensor, indices: Tensor, axis: int) -> Tensor:
@@ -383,7 +364,6 @@ def gather(self, array: Tensor, indices: Tensor, axis: int) -> Tensor:
         Returns:
             array: values of the array at the given indices
         """
-        ...
 
     @abstractmethod
     def hash_tensor(self, tensor: Tensor) -> int:
@@ -395,7 +375,6 @@ def hash_tensor(self, tensor: Tensor) -> int:
         Returns:
             int: hash of the given tensor
         """
-        ...
 
     @abstractmethod
     def hermite_renormalized(self, A: Matrix, B: Vector, C: Scalar, shape: Sequence[int]) -> Tensor:
@@ -410,7 +389,6 @@ def hermite_renormalized(self, A: Matrix, B: Vector, C: Scalar, shape: Sequence[
         Returns:
             array: renormalized hermite polynomials
         """
-        ...
 
     @abstractmethod
     def imag(self, array: Tensor) -> Tensor:
@@ -422,7 +400,6 @@ def imag(self, array: Tensor) -> Tensor:
         Returns:
             array: imaginary part of array
         """
-        ...
 
     @abstractmethod
     def inv(self, tensor: Tensor) -> Tensor:
@@ -434,12 +411,10 @@ def inv(self, tensor: Tensor) -> Tensor:
         Returns:
             array: inverse of tensor
         """
-        ...
 
     @abstractmethod
     def is_trainable(self, tensor: Tensor) -> bool:
         r"""Returns whether the given tensor is trainable."""
-        ...
 
     @abstractmethod
     def lgamma(self, x: Tensor) -> Tensor:
@@ -451,7 +426,6 @@ def lgamma(self, x: Tensor) -> Tensor:
         Returns:
             array: natural logarithm of the gamma function of ``x``
         """
-        ...
 
     @abstractmethod
     def log(self, x: Tensor) -> Tensor:
@@ -463,7 +437,6 @@ def log(self, x: Tensor) -> Tensor:
         Returns:
             array: natural logarithm of ``x``
         """
-        ...
 
     @abstractmethod
     def matmul(
@@ -488,7 +461,6 @@ def matmul(
         Returns:
             array: matrix product of ``a`` and ``b``
         """
-        ...
 
     @abstractmethod
     def matvec(self, a: Matrix, b: Vector, transpose_a=False, adjoint_a=False) -> Tensor:
@@ -503,7 +475,6 @@ def matvec(self, a: Matrix, b: Vector, transpose_a=False, adjoint_a=False) -> Te
         Returns:
             array: matrix vector product of ``a`` and ``b``
         """
-        ...
 
     @abstractmethod
     def maximum(self, a: Tensor, b: Tensor) -> Tensor:
@@ -516,7 +487,6 @@ def maximum(self, a: Tensor, b: Tensor) -> Tensor:
         Returns:
             array: element-wise maximum of ``a`` and ``b``
         """
-        ...
 
     @abstractmethod
     def minimum(self, a: Tensor, b: Tensor) -> Tensor:
@@ -529,7 +499,6 @@ def minimum(self, a: Tensor, b: Tensor) -> Tensor:
         Returns:
             array: element-wise minimum of ``a`` and ``b``
         """
-        ...
 
     @abstractmethod
     def new_variable(
@@ -545,7 +514,6 @@ def new_variable(
         Returns:
             array: new variable
         """
-        ...
 
     @abstractmethod
     def new_constant(self, value: Tensor, name: str, dtype: Any) -> Tensor:
@@ -559,7 +527,6 @@ def new_constant(self, value: Tensor, name: str, dtype: Any) -> Tensor:
         Returns:
             array: new constant
         """
-        ...
 
     @abstractmethod
     def norm(self, array: Tensor) -> Tensor:
@@ -571,7 +538,6 @@ def norm(self, array: Tensor) -> Tensor:
         Returns:
             array: norm of array
         """
-        ...
 
     @abstractmethod
     def ones(self, shape: Sequence[int], dtype) -> Tensor:
@@ -584,7 +550,7 @@ def ones(self, shape: Sequence[int], dtype) -> Tensor:
         Returns:
             array: array of ones
         """
-        ...  # NOTE : should be float64 by default
+        # NOTE : should be float64 by default
 
     @abstractmethod
     def ones_like(self, array: Tensor) -> Tensor:
@@ -596,7 +562,6 @@ def ones_like(self, array: Tensor) -> Tensor:
         Returns:
             array: array of ones
         """
-        ...
 
     @abstractmethod
     def outer(self, array1: Tensor, array2: Tensor) -> Tensor:
@@ -609,7 +574,6 @@ def outer(self, array1: Tensor, array2: Tensor) -> Tensor:
         Returns:
             array: outer product of array1 and array2
         """
-        ...
 
     @abstractmethod
     def pad(
@@ -626,7 +590,6 @@ def pad(
         Returns:
             array: padded array
         """
-        ...
 
     @abstractmethod
     def pinv(self, matrix: Tensor) -> Tensor:
@@ -638,7 +601,6 @@ def pinv(self, matrix: Tensor) -> Tensor:
         Returns:
             array: pseudo-inverse of matrix
         """
-        ...
 
     @abstractmethod
     def pow(self, x: Tensor, y: Tensor) -> Tensor:
@@ -650,7 +612,6 @@ def pow(self, x: Tensor, y: Tensor) -> Tensor:
         Returns:
             array: :math:`x^y`
         """
-        ...
 
     @abstractmethod
     def real(self, array: Tensor) -> Tensor:
@@ -662,7 +623,6 @@ def real(self, array: Tensor) -> Tensor:
         Returns:
             array: real part of ``array``
         """
-        ...
 
     @abstractmethod
     def reshape(self, array: Tensor, shape: Sequence[int]) -> Tensor:
@@ -675,7 +635,6 @@ def reshape(self, array: Tensor, shape: Sequence[int]) -> Tensor:
         Returns:
             array: reshaped array
         """
-        ...
 
     @abstractmethod
     def sin(self, array: Tensor) -> Tensor:
@@ -687,7 +646,6 @@ def sin(self, array: Tensor) -> Tensor:
         Returns:
             array: sine of ``array``
         """
-        ...
 
     @abstractmethod
     def sinh(self, array: Tensor) -> Tensor:
@@ -699,7 +657,6 @@ def sinh(self, array: Tensor) -> Tensor:
         Returns:
             array: hyperbolic sine of ``array``
         """
-        ...
 
     @abstractmethod
     def sqrt(self, x: Tensor, dtype=None) -> Tensor:
@@ -712,7 +669,10 @@ def sqrt(self, x: Tensor, dtype=None) -> Tensor:
         Returns:
             array: square root of ``x``
         """
-        ...
+
+    @abstractmethod
+    def sqrtm(self, tensor: Tensor) -> Tensor:
+        r"""Returns the matrix square root."""
 
     @abstractmethod
     def sum(self, array: Tensor, axes: Sequence[int] = None):
@@ -725,7 +685,6 @@ def sum(self, array: Tensor, axes: Sequence[int] = None):
         Returns:
             array: sum of array
         """
-        ...
 
     @abstractmethod
     def tensordot(self, a: Tensor, b: Tensor, axes: Sequence[int]) -> Tensor:
@@ -739,7 +698,6 @@ def tensordot(self, a: Tensor, b: Tensor, axes: Sequence[int]) -> Tensor:
         Returns:
             array: tensordot product of ``a`` and ``b``
         """
-        ...
 
     @abstractmethod
     def tile(self, array: Tensor, repeats: Sequence[int]) -> Tensor:
@@ -752,7 +710,6 @@ def tile(self, array: Tensor, repeats: Sequence[int]) -> Tensor:
         Returns:
             array: tiled array
         """
-        ...
 
     @abstractmethod
     def trace(self, array: Tensor, dtype: Any = None) -> Tensor:
@@ -765,7 +722,6 @@ def trace(self, array: Tensor, dtype: Any = None) -> Tensor:
         Returns:
             array: trace of array
         """
-        ...
 
     @abstractmethod
     def transpose(self, a: Tensor, perm: Sequence[int] = None):
@@ -778,7 +734,6 @@ def transpose(self, a: Tensor, perm: Sequence[int] = None):
         Returns:
             array: transposed array
         """
-        ...
 
     @abstractmethod
     def unique_tensors(self, lst: List[Tensor]) -> List[Tensor]:
@@ -790,7 +745,6 @@ def unique_tensors(self, lst: List[Tensor]) -> List[Tensor]:
         Returns:
             list: list of tensors without duplicates and non-tensors.
         """
-        ...
 
     @abstractmethod
     def update_tensor(self, tensor: Tensor, indices: Tensor, values: Tensor) -> Tensor:
@@ -801,7 +755,6 @@ def update_tensor(self, tensor: Tensor, indices: Tensor, values: Tensor) -> Tens
             indices (array): indices to update
             values (array): values to update
         """
-        ...
 
     @abstractmethod
     def update_add_tensor(self, tensor: Tensor, indices: Tensor, values: Tensor) -> Tensor:
@@ -812,7 +765,6 @@ def update_add_tensor(self, tensor: Tensor, indices: Tensor, values: Tensor) ->
             indices (array): indices to update
             values (array): values to add
         """
-        ...
 
     @abstractmethod
     def value_and_gradients(
@@ -827,7 +779,6 @@ def value_and_gradients(
         Returns:
             tuple: loss and gradients (dict) of the given cost function
         """
-        ...
 
     @abstractmethod
     def zeros(self, shape: Sequence[int], dtype) -> Tensor:
@@ -840,7 +791,6 @@ def zeros(self, shape: Sequence[int], dtype) -> Tensor:
         Returns:
             array: array of zeros
         """
-        ...
 
     @abstractmethod
     def zeros_like(self, array: Tensor) -> Tensor:
@@ -852,7 +802,6 @@ def zeros_like(self, array: Tensor) -> Tensor:
         Returns:
             array: array of zeros
         """
-        ...
 
     # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     # Methods that build on the basic ops and don't need to be overridden in the backend implementation
@@ -865,16 +814,6 @@ def euclidean_opt(self):
             self._euclidean_opt = self.DefaultEuclideanOptimizer()
         return self._euclidean_opt
 
-    # pylint: disable=no-self-use
-    def eigvals(self, tensor: Tensor) -> Tensor:
-        r"""Returns the eigenvalues of a matrix."""
-        ...
-
-    # pylint: disable=no-self-use
-    def sqrtm(self, tensor: Tensor) -> Tensor:
-        r"""Returns the matrix square root."""
-        ...
-
     # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     # Methods that build on the basic ops and don't need to be overridden in the backend implementation
     # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -916,7 +855,7 @@ def unitary_to_orthogonal(self, U):
         Y = self.imag(U)
         return self.block([[X, -Y], [Y, X]])
 
-    def random_symplectic(self, num_modes: int = 1, max_r: float = 1.0) -> Tensor:
+    def random_symplectic(self, num_modes: int, max_r: float = 1.0) -> Tensor:
         r"""A random symplectic matrix in ``Sp(2*num_modes)``.
 
         Squeezing is sampled uniformly from 0.0 to ``max_r`` (1.0 by default).
@@ -933,13 +872,18 @@ def random_symplectic(self, num_modes: int = 1, max_r: float = 1.0) -> Tensor:
         dd = self.diag(self.concat([self.exp(-r), np.exp(r)], axis=0), k=0)
         return OW @ dd @ OV
 
-    def random_orthogonal(self, num_modes: int = 1) -> Tensor:
-        """A random orthogonal matrix in :math:`O(2*num_modes)`."""
-        if num_modes == 1:
-            W = self.exp(1j * np.random.uniform(size=(1, 1)))
-        else:
-            W = unitary_group.rvs(dim=num_modes)
-        return self.unitary_to_orthogonal(W)
+    @staticmethod
+    def random_orthogonal(N: int) -> Tensor:
+        """A random orthogonal matrix in :math:`O(N)`."""
+        if N == 1:
+            return np.array([[1.0]])
+        return ortho_group.rvs(dim=N)
+
+    def random_unitary(self, N: int) -> Tensor:
+        """a random unitary matrix in :math:`U(N)`"""
+        if N == 1:
+            return self.exp(1j * np.random.uniform(size=(1, 1)))
+        return unitary_group.rvs(dim=N)
 
     def single_mode_to_multimode_vec(self, vec, num_modes: int):
         r"""Apply the same 2-vector (i.e. single-mode) to a larger number of modes."""
diff --git a/mrmustard/utils/training.py b/mrmustard/utils/training.py
index 55a895dcb..ffff5e993 100644
--- a/mrmustard/utils/training.py
+++ b/mrmustard/utils/training.py
@@ -119,55 +119,6 @@ def _render_traceback_(self):
 # ~~~~~~~~~~~~~~~~~
 
 
-# def new_variable(
-#     value, bounds: Tuple[Optional[float], Optional[float]], name: str, dtype=math.float64
-# ) -> Trainable:
-#     r"""Returns a new trainable variable from the current math backend
-#     with initial value set by `value` and bounds set by `bounds`.
-#
-#     Args:
-#         value (float): The initial value of the variable
-#         bounds (Tuple[float, float]): The bounds of the variable
-#         name (str): The name of the variable
-#         dtype: The dtype of the variable
-#
-#     Returns:
-#         variable (Trainable): The new variable
-#     """
-#     return math.new_variable(value, bounds, name, dtype)
-
-
-# def new_constant(value, name: str, dtype=math.float64) -> Tensor:
-#     r"""Returns a new constant (non-trainable) tensor from the current math backend
-#     with initial value set by `value`.
-#     Args:
-#         value (numeric): The initial value of the tensor
-#         name (str): The name of the constant
-#         dtype: The dtype of the constant
-#
-#     Returns:
-#         tensor (Tensor): The new constant tensor
-#     """
-#     return math.new_constant(value, name, dtype)
-
-
-def new_symplectic(num_modes: int) -> Tensor:
-    r"""Returns a new symplectic matrix from the current math backend with ``num_modes`` modes.
-
-    Args:
-        num_modes (int): the number of modes in the symplectic matrix
-
-    Returns:
-        Tensor: the new symplectic matrix
-    """
-    return math.random_symplectic(num_modes)
-
-
-def new_orthogonal(num_modes: int) -> Tensor:
-    """Returns a random orthogonal matrix in :math:`O(2*num_modes)`."""
-    return math.random_orthogonal(num_modes)
-
-
 def update_symplectic(
     symplectic_params: Sequence[Trainable], symplectic_grads: Sequence[Tensor], symplectic_lr: float
 ):
diff --git a/tests/test_utils/test_opt.py b/tests/test_utils/test_opt.py
index 7f38850b7..af3bd994e 100644
--- a/tests/test_utils/test_opt.py
+++ b/tests/test_utils/test_opt.py
@@ -20,7 +20,15 @@
 
 from thewalrus.symplectic import two_mode_squeezing
 
-from mrmustard.lab.gates import Sgate, BSgate, S2gate, Ggate, Interferometer, Ggate
+from mrmustard.lab.gates import (
+    Sgate,
+    BSgate,
+    S2gate,
+    Ggate,
+    Interferometer,
+    Ggate,
+    RealInterferometer,
+)
 from mrmustard.lab.circuit import Circuit
 from mrmustard.utils.training import Optimizer
 from mrmustard.utils.parametrized import Parametrized
@@ -181,6 +189,31 @@ def cost_fn():
     assert np.allclose(-cost_fn(), 0.25, atol=1e-5)
 
 
+def test_learning_two_mode_RealInterferometer():
+    """Finding the optimal Interferometer to make a pair of single photons"""
+    np.random.seed(11)
+    ops = [
+        Sgate(
+            r=np.random.normal(size=(2)) ** 2,
+            phi=np.random.normal(size=(2)),
+            r_trainable=True,
+            phi_trainable=True,
+        ),
+        RealInterferometer(num_modes=2, orthogonal_trainable=True),
+    ]
+    circ = Circuit(ops)
+    state_in = Vacuum(num_modes=2)
+
+    def cost_fn():
+        amps = (state_in >> circ).ket(cutoffs=[2, 2])
+        return -tf.abs(amps[1, 1]) ** 2 + tf.abs(amps[0, 1]) ** 2
+
+    opt = Optimizer(orthogonal_lr=0.5, euclidean_lr=0.01)
+
+    opt.minimize(cost_fn, by_optimizing=[circ], max_steps=1000)
+    assert np.allclose(-cost_fn(), 0.25, atol=1e-5)
+
+
 def test_learning_four_mode_Interferometer():
     """Finding the optimal Interferometer to make a NOON state with N=2"""
     np.random.seed(11)
@@ -214,6 +247,39 @@ def cost_fn():
     assert np.allclose(-cost_fn(), 0.0625, atol=1e-5)
 
 
+def test_learning_four_mode_RealInterferometer():
+    """Finding the optimal Interferometer to make a NOON state with N=2"""
+    np.random.seed(11)
+    ops = [
+        Sgate(
+            r=np.random.uniform(size=4),
+            phi=np.random.normal(size=4),
+            r_trainable=True,
+            phi_trainable=True,
+        ),
+        RealInterferometer(num_modes=4, orthogonal_trainable=True),
+    ]
+    circ = Circuit(ops)
+    state_in = Vacuum(num_modes=4)
+
+    def cost_fn():
+        amps = (state_in >> circ).ket(cutoffs=[3, 3, 3, 3])
+        return (
+            -tf.abs(
+                tf.reduce_sum(
+                    amps[1, 1]
+                    * np.array([[0, 0, 1 / np.sqrt(2)], [0, 0, 0], [1 / np.sqrt(2), 0, 0]])
+                )
+            )
+            ** 2
+        )
+
+    opt = Optimizer(symplectic_lr=0.5, euclidean_lr=0.01)
+
+    opt.minimize(cost_fn, by_optimizing=[circ], max_steps=1000)
+    assert np.allclose(-cost_fn(), 0.0625, atol=1e-5)
+
+
 def test_squeezing_hong_ou_mandel_optimizer():
     """Finding the optimal squeezing parameter to get Hong-Ou-Mandel dip in time
     see https://www.pnas.org/content/117/52/33107/tab-article-info