Implements the Montrealer

.github/ACKNOWLEDGMENTS.md

@@ -67,3 +67,5 @@
 * [Martin Houde](https://github.com/MHoude2) (Polytechnique Montréal) - 🙃 Minister of amplification
 
 * Will McCutcheon (Heriot-Watt University) - 🧅 Gaussian Onion Merchant
+
+* [Yanic Cardin](https://github.com/yaniccd) (Polytechnique Montréal) - 🦜 Pirate of the permutations

.github/CHANGELOG.md

@@ -2,8 +2,9 @@
 
 ### New features
 
-* Adds the Takagi decomposition [(#363)](https://github.com/XanaduAI/thewalrus/pull/338)
+* Adds the Takagi decomposition [(#363)](https://github.com/XanaduAI/thewalrus/pull/363)
 
+* Adds the Montrealer and Loop Montrealer functions [(#363)](https://github.com/XanaduAI/thewalrus/pull/374).
 ### Breaking changes
 
 ### Improvements
@@ -25,7 +26,7 @@
 
 This release contains contributions from (in alphabetical order): 
 
-Gregory Morse, Nicolas Quesada
+Yanic Cardin, Gregory Morse, Nicolas Quesada
 
 ---
 

thewalrus/__init__.py

@@ -132,6 +132,12 @@
     rec_torontonian,
     rec_ltorontonian,
 )
+
+from ._montrealer import (
+    mtl,
+    lmtl,
+)
+
 from ._version import __version__
 
 
@@ -152,6 +158,8 @@
     "reduction",
     "hermite_multidimensional",
     "grad_hermite_multidimensional",
+    "mtl",
+    "lmtl",
     "version",
 ]
 

thewalrus/_montrealer.py

@@ -0,0 +1,135 @@
+"""
+Montrealer Python interface
+
+* Yanic Cardin and Nicolás Quesada. "Photon-number moments and cumulants of Gaussian states"
+  `arxiv:12212.06067 (2023) <https://arxiv.org/abs/arxiv:2212.06067v2>`_
+"""
+import numpy as np
+import numba
+from thewalrus.quantum.conversions import Xmat
+from thewalrus.charpoly import powertrace
+from ._hafnian import nb_ix
+from ._torontonian import tor_input_checks
+
+
+@numba.jit(nopython=True, cache=True)
+def dec2bin(num, n):  # pragma: no cover
+    """Helper function to generate convert an integer into an element of the powerset of n objects.
+
+    Args:
+        num (int): label to convert.
+        n (int): number of elements in the set.
+
+    Returns:
+        (array): array containing the labels of the elements to be selected.
+
+    """
+    digits = np.zeros((n), dtype=type(num))
+    nn = num
+    counter = -1
+    while nn >= 1:  # Should be >= 1, not > 1.
+        digits[counter] = nn % 2
+        counter -= 1
+        nn //= 2
+    return np.nonzero(digits)[0]
+
+
+@numba.jit(nopython=True)
+def montrealer(Sigma):  # pragma: no cover
+    """Calculates the loop-montrealer of the zero-displacement Gaussian state with the given complex covariance matrix.
+
+    Args:
+        A (array): adjacency matrix of the Gaussian state.
+
+    Returns:
+        (np.complex128): the montrealer of A.
+
+    """
+    n = len(Sigma) // 2
+    tot_num = 2**n
+    val = np.complex128(0)
+    for p in numba.prange(tot_num):
+        pos = dec2bin(p, n)
+        lenpos = len(pos)
+        pos = np.concatenate((pos, n + pos))
+        submat = nb_ix(Sigma, pos, pos)
+        sign = (-1) ** (lenpos + 1)
+        val += (sign) * powertrace(submat, n + 1)[-1]
+    return (-1) ** (n + 1) * val / (2 * n)
+
+
+@numba.jit(nopython=True)
+def power_loop(Sigma, zeta, n):  # pragma: no cover
+    """Auxiliary function to calculate the product np.conj(zeta) @ Sigma^{n-1} @ zeta.
+    Args:
+        Sigma (array): square complex matrix
+        zeta (array): complex vector
+        n (int): sought after power
+
+    Returns:
+        (np.complex128 or np.float64): the product np.conj(zeta) @ Sigma^{n-1} @ zeta
+    """
+    vec = zeta
+    for _ in range(n - 1):
+        vec = Sigma @ vec
+    return np.conj(zeta) @ vec
+
+
+@numba.jit(nopython=True, cache=True)
+def lmontrealer(Sigma, zeta):  # pragma: no cover
+    """Calculates the loop-montrealer of the displaced Gaussian state with the given complex covariance matrix and vector of displacements.
+
+    Args:
+        Sigma (array): complex Glauber covariance matrix of the Gaussian state.
+        zeta (array): vector of displacements
+
+    Returns:
+        (np.complex128): the montrealer of A.
+    """
+    n = len(Sigma) // 2
+    tot_num = 2**n
+    val = np.complex128(0)
+    val_loops = np.complex128(0)
+    for p in numba.prange(tot_num):
+        pos = dec2bin(p, n)
+        lenpos = len(pos)
+        pos = np.concatenate((pos, n + pos))
+        subvec = zeta[pos]
+        submat = nb_ix(Sigma, pos, pos)
+        sign = (-1) ** (lenpos + 1)
+        val_loops += sign * power_loop(submat, subvec, n)
+        val += sign * powertrace(submat, n + 1)[-1]
+    return (-1) ** (n + 1) * (val / (2 * n) + val_loops / 2)
+
+
+def lmtl(A, zeta):
+    """Returns the montrealer of an NxN matrix and an N-length vector.
+
+    Args:
+        A (array): an NxN array of even dimensions.
+        zeta (array): an N-length vector of even dimensions
+
+    Returns:
+        np.float64 or np.complex128: the loop montrealer of matrix A, vector zeta
+    """
+
+    tor_input_checks(A, zeta)
+    n = len(A) // 2
+    Sigma = Xmat(n) @ A
+    return lmontrealer(Sigma, zeta)
+
+
+def mtl(A):
+    """Returns the montrealer of an NxN matrix.
+
+    Args:
+        A (array): an NxN array of even dimensions.
+
+    Returns:
+        np.float64 or np.complex128: the montrealer of matrix A
+    """
+
+    tor_input_checks(A)
+    n = len(A) // 2
+    Sigma = Xmat(n) @ A
+    return montrealer(Sigma)

thewalrus/_torontonian.py

@@ -20,26 +20,40 @@
 from ._hafnian import reduction, find_kept_edges, nb_ix
 
 
-def tor(A, recursive=True):
-    """Returns the Torontonian of a matrix.
-
+def tor_input_checks(A, loops=None):
+    """Checks the correcteness of the inputs for the torontonian/montrealer.
     Args:
-        A (array): a square array of even dimensions.
-        recursive: use the faster recursive implementation.
-
-    Returns:
-        np.float64 or np.complex128: the torontonian of matrix A.
+        A (array): an NxN array of even dimensions.
+        loops (array): optinal argument, an N-length vector of even dimensions.
     """
     if not isinstance(A, np.ndarray):
         raise TypeError("Input matrix must be a NumPy array.")
-
     matshape = A.shape
 
     if matshape[0] != matshape[1]:
         raise ValueError("Input matrix must be square.")
 
     if matshape[0] % 2 != 0:
         raise ValueError("matrix dimension must be even")
+
+    if loops is not None:
+        if not isinstance(loops, np.ndarray):
+            raise TypeError("Input matrix must be a NumPy array.")
+        if matshape[0] != len(loops):
+            raise ValueError("gamma must be a vector matching the dimension of A")
+
+
+def tor(A, recursive=True):
+    """Returns the Torontonian of a matrix.
+
+    Args:
+        A (array): a square array of even dimensions.
+        recursive: use the faster recursive implementation.
+
+    Returns:
+        np.float64 or np.complex128: the torontonian of matrix A.
+    """
+    tor_input_checks(A)
     return rec_torontonian(A) if recursive else numba_tor(A)
 
 
@@ -54,23 +68,7 @@ def ltor(A, gamma, recursive=True):
     Returns:
         np.float64 or np.complex128: the loop torontonian of matrix A, vector gamma
     """
-
-    if not isinstance(A, np.ndarray):
-        raise TypeError("Input matrix must be a NumPy array.")
-
-    if not isinstance(gamma, np.ndarray):
-        raise TypeError("Input matrix must be a NumPy array.")
-
-    matshape = A.shape
-
-    if matshape[0] != matshape[1]:
-        raise ValueError("Input matrix must be square.")
-
-    if matshape[0] != len(gamma):
-        raise ValueError("gamma must be a vector matching the dimension of A")
-
-    if matshape[0] % 2 != 0:
-        raise ValueError("matrix dimension must be even")
+    tor_input_checks(A, gamma)
 
     return rec_ltorontonian(A, gamma) if recursive else numba_ltor(A, gamma)