feat(gaussian): Loop torontonian implementation

The loop torontonian is implemented according to
https://arxiv.org/abs/2202.04600 and https://arxiv.org/abs/2109.04528.

Moreover, the OMP race condition is fixed in `torontonian.cpp`.

CMakeLists.txt

@@ -16,9 +16,10 @@ set(CMAKE_CXX_STANDARD 17)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
 set(CMAKE_POSITION_INDEPENDENT_CODE ON)
 
-add_subdirectory(src)
-add_subdirectory(piquasso/_math)
-
 if(CMAKE_BUILD_TYPE STREQUAL "Debug")
     add_definitions(-DDEBUG)
 endif()
+
+
+add_subdirectory(src)
+add_subdirectory(piquasso/_math)

piquasso/_math/CMakeLists.txt

@@ -1,5 +1,5 @@
 pybind11_add_module(torontonian SHARED ${CMAKE_CURRENT_SOURCE_DIR}/torontonian.cpp)
 
-target_link_libraries(torontonian PUBLIC torontonianboost)
+target_link_libraries(torontonian PUBLIC torontonianboost looptorontonianboost)
 
 install(TARGETS torontonian LIBRARY DESTINATION piquasso/_math)

piquasso/_math/linalg.py

@@ -81,6 +81,10 @@ def block_reduce(array: np.ndarray, reduce_on: Tuple[int, ...]) -> np.ndarray:
     return reduce_(array, reduce_on=(reduce_on * 2))
 
 
+def block_reduce_xpxp(array: np.ndarray, reduce_on: Tuple[int, ...]) -> np.ndarray:
+    return reduce_(array, reduce_on=np.repeat(reduce_on, 2))
+
+
 def assym_reduce(array, row_reduce_on, col_reduce_on):
     particles = np.sum(row_reduce_on)
 

piquasso/_math/torontonian.cpp

@@ -19,21 +19,57 @@
 
 #include "numpy_utils.hpp"
 #include "torontonian.hpp"
+#include "loop_torontonian.hpp"
 
 namespace py = pybind11;
 
 template <typename TScalar>
-py::object torontonian_np(py::array_t<TScalar, py::array::c_style | py::array::forcecast> array)
+py::object torontonian_np(
+    py::array_t<TScalar, py::array::c_style | py::array::forcecast> matrix)
 {
-    Matrix<TScalar> matrix = numpy_to_matrix(array);
+    Matrix<TScalar> native_matrix = numpy_to_matrix(matrix);
 
-    TScalar result = torontonian_cpp(matrix);
+    TScalar result = torontonian_cpp(native_matrix);
 
     return create_numpy_scalar(result);
 }
 
+template <typename TScalar>
+py::object loop_torontonian_np(
+    py::array_t<TScalar, py::array::c_style | py::array::forcecast> matrix,
+    py::array_t<TScalar, py::array::c_style | py::array::forcecast> displacement_vector)
+{
+    Matrix<TScalar> native_matrix = numpy_to_matrix(matrix);
+    Vector<TScalar> native_displacement_vector = numpy_to_vector(displacement_vector);
+
+    TScalar result = loop_torontonian_cpp(native_matrix, native_displacement_vector);
+
+    return create_numpy_scalar(result);
+}
+
+const char* torontonian_docstring = R""""(
+Calculates the torontonian of a matrix.
+
+Note:
+    This function expects arguments in a different ordering than usual. Usually, the
+    inputs are in the xxpp-ordering, but for this implementation, one needs to provide
+    data in the xpxp-ordering.
+)"""";
+
+
+const char* loop_torontonian_docstring = R""""(
+Calculates the loop torontonian of a matrix and a displacement vector.
+
+Note:
+    This function expects arguments in a different ordering than usual. Usually, the
+    inputs are in the xxpp-ordering, but for this implementation, one needs to provide
+    data in the xpxp-ordering.
+)"""";
+
 PYBIND11_MODULE(torontonian, m)
 {
-    m.def("torontonian", &torontonian_np<float>);
-    m.def("torontonian", &torontonian_np<double>);
+    m.def("torontonian", &torontonian_np<float>, torontonian_docstring);
+    m.def("torontonian", &torontonian_np<double>, torontonian_docstring);
+    m.def("loop_torontonian", &loop_torontonian_np<float>, loop_torontonian_docstring);
+    m.def("loop_torontonian", &loop_torontonian_np<double>, loop_torontonian_docstring);
 }

piquasso/_simulators/gaussian/calculations.py

@@ -25,7 +25,10 @@
 from scipy.special import factorial
 
 from .state import GaussianState
-from .probabilities import calculate_click_probability
+from .probabilities import (
+    calculate_click_probability_nondisplaced,
+    calculate_click_probability,
+)
 
 from piquasso.instructions import gates
 
@@ -512,18 +515,11 @@ def threshold_measurement(
     return Result(state=state, samples=samples)
 
 
-def _generate_threshold_samples_using_torontonian(
-    state,
-    instruction,
-    shots,
-):
-    if not np.allclose(state.xpxp_mean_vector, np.zeros_like(state.xpxp_mean_vector)):
-        raise NotImplementedError(
-            "Threshold measurement for displaced states are not supported: "
-            f"xpxp_mean_vector={state.xpxp_mean_vector}"
-        )
+def _generate_threshold_samples_using_torontonian(state, instruction, shots):
+    is_displaced = state._is_displaced()
 
     rng = state._config.rng
+    hbar = state._config.hbar
 
     modes = instruction.modes
 
@@ -533,10 +529,16 @@ def get_probability(
     ) -> float:
         reduced_state = state.reduced(subspace_modes)
 
+        if not is_displaced:
+            return calculate_click_probability_nondisplaced(
+                reduced_state.xpxp_covariance_matrix / hbar,
+                occupation_numbers,
+            )
+
         return calculate_click_probability(
-            reduced_state.xxpp_covariance_matrix,
+            reduced_state.xpxp_covariance_matrix / hbar,
+            reduced_state.xpxp_mean_vector / np.sqrt(hbar),
             occupation_numbers,
-            state._config.hbar,
         )
 
     samples = []

piquasso/_simulators/gaussian/probabilities.py

@@ -21,8 +21,8 @@
 
 from scipy.special import factorial
 
-from piquasso._math.linalg import block_reduce
-from piquasso._math.torontonian import torontonian
+from piquasso._math.linalg import block_reduce_xpxp
+from piquasso._math.torontonian import torontonian, loop_torontonian
 from piquasso.api.calculator import BaseCalculator
 
 
@@ -142,10 +142,9 @@ def calculate_hafnian(self, reduce_on: np.ndarray) -> float:
         return self.calculator.loop_hafnian(self._A, self._gamma, reduce_on)
 
 
-def calculate_click_probability(
-    xp_covariance: np.ndarray,
+def calculate_click_probability_nondisplaced(
+    xpxp_covariance: np.ndarray,
     occupation_number: Tuple[int, ...],
-    hbar: float,
 ) -> float:
     r"""
     Calculates the threshold detection probability with the equation
@@ -167,19 +166,71 @@ def calculate_click_probability(
             \right ).
     """
 
-    d = len(xp_covariance) // 2
+    d = len(xpxp_covariance) // 2
 
-    sigma: np.ndarray = (xp_covariance / hbar + np.identity(2 * d)) / 2
+    sigma: np.ndarray = (xpxp_covariance + np.identity(2 * d)) / 2
 
-    sigma_inv_reduced = block_reduce(
-        np.linalg.inv(sigma),
-        reduce_on=occupation_number,
-    )
+    sigma_inv = np.linalg.inv(sigma)
+
+    sigma_inv_reduced = block_reduce_xpxp(sigma_inv, reduce_on=occupation_number)
+
+    A = np.identity(len(sigma_inv_reduced), dtype=float) - sigma_inv_reduced
+
+    probability = torontonian(A) / np.sqrt(np.linalg.det(sigma))
+
+    return max(probability, 0.0)
+
+
+def calculate_click_probability(
+    xpxp_covariance: np.ndarray,
+    xpxp_mean: np.ndarray,
+    occupation_number: Tuple[int, ...],
+) -> float:
+    r"""
+    Calculates the threshold detection probability with the equation
+
+    .. math::
+        p(S) = \frac{
+            \operatorname{ltor}( I - ( \Sigma^{-1} )_{(S)}, \vec{\gamma} )
+        }{
+            \sqrt{\operatorname{det}(\Sigma)}
+        },
+
+    where :math:`\Sigma \in \mathbb{R}^{2d \times 2d}` is a symmetric matrix
+    defined by
+
+    .. math::
+        \Sigma = \frac{1}{2} \left (
+                \frac{1}{\hbar} \sigma_{xp}
+                + I
+            \right ),
+
+    and
+
+    .. math::
+        \vec{\gamma} = \Sigma^{-1} \vec{\alpha}
+    """
+
+    d = len(xpxp_covariance) // 2
+
+    sigma: np.ndarray = (xpxp_covariance + np.identity(2 * d)) / 2
+
+    sigma_inv = np.linalg.inv(sigma)
+
+    gamma = sigma_inv @ xpxp_mean
+
+    sigma_inv_reduced = block_reduce_xpxp(sigma_inv, reduce_on=occupation_number)
+
+    gamma_reduced = block_reduce_xpxp(gamma, reduce_on=occupation_number)
+
+    A = np.identity(len(sigma_inv_reduced), dtype=float) - sigma_inv_reduced
+
+    exponential_term = np.exp(-xpxp_mean @ sigma_inv @ xpxp_mean / 2)
 
     probability = (
-        torontonian(
-            np.identity(len(sigma_inv_reduced), dtype=float) - sigma_inv_reduced
-        )
-    ).real / np.sqrt(np.linalg.det(sigma))
+        loop_torontonian(A, gamma_reduced).real
+        * exponential_term
+        / np.sqrt(np.linalg.det(sigma))
+    )
 
     return max(probability, 0.0)

piquasso/_simulators/gaussian/state.py

@@ -39,6 +39,7 @@
     DensityMatrixCalculation,
     DisplacedDensityMatrixCalculation,
     NondisplacedDensityMatrixCalculation,
+    calculate_click_probability_nondisplaced,
     calculate_click_probability,
 )
 
@@ -929,10 +930,18 @@ def get_threshold_detection_probability(
                 f"occupation_number='{occupation_number}'."
             )
 
+        hbar = self._config.hbar
+
+        if not self._is_displaced():
+            return calculate_click_probability_nondisplaced(
+                self.xpxp_covariance_matrix / hbar,
+                tuple(occupation_number),
+            )
+
         return calculate_click_probability(
-            self.xxpp_covariance_matrix,
+            self.xpxp_covariance_matrix / hbar,
+            self.xpxp_mean_vector / np.sqrt(hbar),
             tuple(occupation_number),
-            self._config.hbar,
         )
 
     @property

scripts/loop_torontonian_benchmark.py

@@ -0,0 +1,99 @@
+#
+# Copyright 2021-2024 Budapest Quantum Computing Group
+#
+# 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.
+
+"""
+Benchmarking the Piquasso and TheWalrus torontonian implementations.
+"""
+
+import time
+
+import piquasso as pq
+
+from piquasso._math.torontonian import loop_torontonian as pq_loop_torontonian
+from thewalrus import ltor as tw_loop_torontonian
+
+import matplotlib.pyplot as plt
+
+from scipy.stats import unitary_group
+
+import numpy as np
+
+np.set_printoptions(suppress=True, linewidth=200)
+
+
+if __name__ == "__main__":
+    x = []
+    y = []
+    z = []
+
+    ITER = 10
+
+    for d in range(3, 20):
+        print(d)
+        x.append(d)
+
+        simulator = pq.GaussianSimulator(d=d)
+
+        program = pq.Program(
+            instructions=[pq.Vacuum()]
+            + [pq.Displacement(r=np.random.rand()).on_modes(i) for i in range(d)]
+            + [pq.Squeezing(r=np.random.rand()).on_modes(i) for i in range(d)]
+            + [pq.Interferometer(unitary_group.rvs(d))]
+        )
+
+        state = simulator.execute(program).state
+
+        xpxp_covariance_matrix = state.xpxp_covariance_matrix
+
+        sigma: np.ndarray = (xpxp_covariance_matrix / 2 + np.identity(2 * d)) / 2
+
+        input_matrix = np.identity(len(sigma), dtype=float) - np.linalg.inv(sigma)
+        displacement_vector = state.xpxp_mean_vector
+
+        sum_ = 0.0
+
+        for _ in range(ITER):
+            print("|", end="", flush=True)
+            start_time = time.time()
+            pq_loop_torontonian(input_matrix, displacement_vector)
+            sum_ += time.time() - start_time
+
+        y.append(sum_ / ITER)
+        print()
+
+        xxpp_covariance_matrix = state.xxpp_covariance_matrix
+
+        sigma: np.ndarray = (xxpp_covariance_matrix / 2 + np.identity(2 * d)) / 2
+
+        input_matrix = np.identity(len(sigma), dtype=float) - np.linalg.inv(sigma)
+        displacement_vector = state.xxpp_mean_vector
+
+        sum_ = 0.0
+
+        for _ in range(ITER):
+            print("-", end="", flush=True)
+            start_time = time.time()
+            tw_loop_torontonian(input_matrix, displacement_vector)
+            sum_ += time.time() - start_time
+
+        z.append(sum_ / ITER)
+        print()
+
+    plt.scatter(x[1:], y[1:], label="Piquasso")
+    plt.scatter(x[1:], z[1:], label="TheWalrus")
+    plt.legend()
+    plt.yscale("log")
+
+    plt.show()