perf(fock): Euler decomposition

`FockSpace.get_linear_fock_operator` has been refactored to avoid embedding
into Fock space. For clarity, the Euler decomposition has been re-implemented
in `euler` in `decompositions`.

piquasso/_backends/fock/general/calculations.py

@@ -22,6 +22,8 @@
 
 from piquasso.instructions import gates
 
+from piquasso._math.decompositions import euler
+
 from piquasso.api.instruction import Instruction
 from piquasso.api.result import Result
 
@@ -41,18 +43,22 @@ def passive_linear(
         state._calculator, state._config
     )
 
+    _apply_passive_linear(state, operator, instruction.modes)
+
+    return Result(state=state)
+
+
+def _apply_passive_linear(state, interferometer, modes):
     fock_operator = state._space.get_passive_fock_operator(
-        operator,
-        modes=instruction.modes,
+        interferometer,
+        modes=modes,
         d=state._space.d,
     )
 
     state._density_matrix = (
         fock_operator @ state._density_matrix @ fock_operator.conjugate().transpose()
     )
 
-    return Result(state=state)
-
 
 def particle_number_measurement(
     state: FockState, instruction: Instruction, shots: int
@@ -290,16 +296,31 @@ def squeezing(state: FockState, instruction: Instruction, shots: int) -> Result:
 def linear(
     state: FockState, instruction: gates._ActiveLinearGate, shots: int
 ) -> Result:
-    operator = state._space.get_linear_fock_operator(
-        modes=instruction.modes,
-        passive_block=instruction._get_passive_block(state._calculator, state._config),
-        active_block=instruction._get_active_block(state._calculator, state._config),
-    )
+    calculator = state._calculator
+    modes = instruction.modes
 
-    state._density_matrix = (
-        operator @ state._density_matrix @ operator.conjugate().transpose()
+    np = calculator.np
+
+    passive_block = instruction._get_passive_block(state._calculator, state._config)
+    active_block = instruction._get_active_block(state._calculator, state._config)
+
+    symplectic = calculator.block(
+        [
+            [passive_block, active_block],
+            [np.conj(active_block), np.conj(passive_block)],
+        ],
     )
 
+    unitary_last, squeezings, unitary_first = euler(symplectic, calculator)
+
+    _apply_passive_linear(state, unitary_first, modes)
+
+    for mode, r in zip(instruction.modes, squeezings):
+        matrix = state._space.get_single_mode_squeezing_operator(r=r, phi=0.0)
+        _apply_active_gate_matrix_to_state(state, matrix, mode)
+
+    _apply_passive_linear(state, unitary_last, modes)
+
     state.normalize()
 
     return Result(state=state)

piquasso/_backends/fock/pure/calculations/__init__.py

@@ -20,11 +20,15 @@
 import random
 import numpy as np
 
+from .passive_linear import _apply_passive_linear
+
 from ...calculations import calculate_state_index_matrix_list
 
 from ..state import PureFockState
 from ..batch_state import BatchPureFockState
 
+from piquasso._math.decompositions import euler
+
 from piquasso.instructions import gates
 
 from piquasso.api.result import Result
@@ -318,13 +322,30 @@ def cubic_phase(state: PureFockState, instruction: Instruction, shots: int) -> R
 def linear(
     state: PureFockState, instruction: gates._ActiveLinearGate, shots: int
 ) -> Result:
-    operator = state._space.get_linear_fock_operator(
-        modes=instruction.modes,
-        passive_block=instruction._get_passive_block(state._calculator, state._config),
-        active_block=instruction._get_active_block(state._calculator, state._config),
+    calculator = state._calculator
+    modes = instruction.modes
+
+    np = calculator.np
+
+    passive_block = instruction._get_passive_block(state._calculator, state._config)
+    active_block = instruction._get_active_block(state._calculator, state._config)
+
+    symplectic = calculator.block(
+        [
+            [passive_block, active_block],
+            [np.conj(active_block), np.conj(passive_block)],
+        ],
     )
 
-    state._state_vector = operator @ state._state_vector
+    unitary_last, squeezings, unitary_first = euler(symplectic, calculator)
+
+    _apply_passive_linear(state, unitary_first, modes, calculator)
+
+    for mode, r in zip(instruction.modes, squeezings):
+        matrix = state._space.get_single_mode_squeezing_operator(r=r, phi=0.0)
+        _apply_active_gate_matrix_to_state(state, matrix, mode)
+
+    _apply_passive_linear(state, unitary_last, modes, calculator)
 
     state.normalize()
 

piquasso/_backends/fock/pure/calculations/passive_linear.py

@@ -44,17 +44,19 @@ def passive_linear(
         state._calculator, state._config
     ).astype(state._config.complex_dtype)
 
+    _apply_passive_linear(state, interferometer, instruction.modes, calculator)
+
+    return Result(state=state)
+
+
+def _apply_passive_linear(state, interferometer, modes, calculator):
     subspace = state._get_subspace(dim=len(interferometer))
 
     subspace_transformations = _get_interferometer_on_fock_space(
         interferometer, subspace, calculator
     )
 
-    _apply_passive_gate_matrix_to_state(
-        state, subspace_transformations, instruction.modes
-    )
-
-    return Result(state=state)
+    _apply_passive_gate_matrix_to_state(state, subspace_transformations, modes)
 
 
 def _get_interferometer_on_fock_space(interferometer, space, calculator):

piquasso/_math/decompositions.py

@@ -84,9 +84,11 @@ def _rotation_to_positive_above_diagonals(block_diagonal_matrix):
 
     return block_diag(
         *[
-            identity
-            if block_diagonal_matrix[2 * index, 2 * index + 1] > 0
-            else rotation
+            (
+                identity
+                if block_diagonal_matrix[2 * index, 2 * index + 1] > 0
+                else rotation
+            )
             for index in range(d)
         ]
     )
@@ -260,3 +262,28 @@ def mean_photon_number_gradient(scaling: float) -> float:
         )
 
     return result.root
+
+
+def euler(symplectic, calculator):
+    np = calculator.np
+    d = len(symplectic) // 2
+
+    identity = np.identity(d)
+    zeros = np.zeros(shape=(d, d), dtype=complex)
+
+    U_orig, R = calculator.polar(symplectic, side="left")
+
+    K = calculator.block(
+        [
+            [identity, zeros],
+            [zeros, -identity],
+        ],
+    )
+
+    H_active = 1j * K @ calculator.logm(R)
+
+    Z = 1j * H_active[:d, d:]
+
+    D, U = takagi(Z, calculator)
+
+    return U, D, np.conj(U).T @ U_orig[:d, :d]

piquasso/_math/fock.py

@@ -24,7 +24,6 @@
 from piquasso.api.config import Config
 from piquasso._math.indices import get_operator_index
 from piquasso._math.combinatorics import partitions
-from piquasso._math.decompositions import takagi
 from piquasso._math.gradients import (
     create_single_mode_displacement_gradient,
     create_single_mode_squeezing_gradient,
@@ -240,148 +239,6 @@ def get_single_mode_cubic_phase_operator(
         position = (annih.T + annih) * np.sqrt(hbar / 2)
         return calculator.expm(1j * calculator.powm(position, 3) * (gamma / (3 * hbar)))
 
-    def embed_matrix(
-        self,
-        matrix: np.ndarray,
-        modes: Tuple[int, ...],
-        auxiliary_modes: Tuple[int, ...],
-    ) -> np.ndarray:
-        indices = []
-        updates = []
-
-        for embedded_index, operator_basis in self.operator_basis_diagonal_on_modes(
-            modes=auxiliary_modes
-        ):
-            index = (
-                operator_basis.ket.on_modes(modes=modes),
-                operator_basis.bra.on_modes(modes=modes),
-            )
-
-            indices.append(embedded_index)
-            updates.append(matrix[index][0])
-
-        return self._calculator.scatter(indices, updates, shape=(self.cardinality,) * 2)
-
-    def get_linear_fock_operator(
-        self,
-        *,
-        modes: Tuple[int, ...],
-        active_block: np.ndarray,
-        passive_block: np.ndarray,
-    ) -> np.ndarray:
-        r"""The matrix of the symplectic transformation in Fock space.
-
-        Any symplectic transformation (in complex representation) can be written as
-
-        .. math::
-            S = \begin{bmatrix}
-                P & A \\
-                A^* & P^*
-            \end{bmatrix}.
-
-        As a first step, this symplectic matrix is polar decomposed:
-
-        .. math::
-            S = R U,
-
-        where :math:`R` is a hermitian matrix and :math:`U` is a unitary one. The polar
-        decomposition of a symplectic matrix is also a symplectic matrix, therefore
-        :math:`U` (being unitary) corresponds to a passive transformation, and :math:`R`
-        to an active one.
-
-        The symplectic matrix :math:`R` has the form
-
-        .. math::
-            \exp \left ( i K \begin{bmatrix}
-                0 & Z \\
-                Z^* & 0
-            \end{bmatrix}
-            \right ),
-
-        where :math:`Z` is a (complex) symmetric matrix. This can be decomposed via
-        Takagi decomposition:
-
-        .. math::
-            Z = U D U^T,
-
-        where :math:`U` is a unitary matrix and :math:`D` is a diagonal matrix. The
-        diagonal entries in :math:`D` correspond to squeezing amplitudes, and :math:`U`
-        corresponds to an interferometer.
-
-        Args:
-            modes (Tuple[int, ...]):
-                The modes on which the transformation should be applied.
-            active_block (np.ndarray): Active part of the symplectic transformation.
-            passive_block (np.ndarray): Passive part of the symplectic transformation.
-
-        Returns:
-            np.ndarray:
-                The resulting transformation, which could be applied to the state.
-        """
-        np = self._calculator.np
-        fallback_np = self._calculator.fallback_np
-
-        d = len(modes)
-        identity = np.identity(d)
-        zeros = np.zeros_like(identity)
-
-        K = self._calculator.block(
-            [
-                [identity, zeros],
-                [zeros, -identity],
-            ],
-        )
-
-        symplectic = self._calculator.block(
-            [
-                [passive_block, active_block],
-                [np.conj(active_block), np.conj(passive_block)],
-            ],
-        )
-
-        U, R = self._calculator.polar(symplectic, side="left")
-
-        H_active = 1j * K @ self._calculator.logm(R)
-        H_passive = 1j * K @ self._calculator.logm(U)
-
-        singular_values, unitary = takagi(1j * H_active[:d, d:], self._calculator)
-
-        transformation = np.identity(self.cardinality, dtype=self.config.complex_dtype)
-
-        fock_operator = self.get_passive_fock_operator(
-            np.conj(unitary).T @ H_passive[:d, :d],
-            modes=modes,
-            d=self.d,
-        )
-
-        transformation = fock_operator @ transformation
-
-        for index, mode in enumerate(modes):
-            operator = self.get_single_mode_squeezing_operator(
-                r=singular_values[index],
-                phi=0.0,
-            )
-
-            squeezing_matrix = self.embed_matrix(
-                operator,
-                modes=(mode,),
-                auxiliary_modes=tuple(
-                    fallback_np.delete(fallback_np.arange(self.d), (mode,))
-                ),
-            )
-
-            transformation = squeezing_matrix @ transformation
-
-        fock_operator = self.get_passive_fock_operator(
-            unitary,
-            modes=modes,
-            d=self.d,
-        )
-
-        transformation = fock_operator @ transformation
-
-        return transformation
-
     @property
     def cardinality(self) -> int:
         return cutoff_cardinality(cutoff=self.cutoff, d=self.d)

tests/backends/test_backend_equivalence.py

@@ -412,7 +412,12 @@ def test_fock_probabilities_with_two_mode_squeezing(SimulatorClass):
 
     probabilities = state.fock_probabilities
 
-    assert all(probability >= 0 for probability in probabilities)
+    # NOTE: It should not happen that the probabilities are negative, but in some cases
+    # (possibly due to floating point errors) they turn out to be slightly negative.
+    assert all(
+        probability >= 0.0 or np.isclose(probability, 0.0)
+        for probability in probabilities
+    )
     assert sum(probabilities) <= 1.0 or np.isclose(sum(probabilities), 1.0)
 
     assert is_proportional(