Hermite improve

docs/quick_guide.rst

@@ -12,7 +12,7 @@ This section provides a quick guide to find which function does what in The Walr
 Numerical hafnian                                                                :func:`thewalrus.hafnian()`
 Symbolic hafnian                                                                 :func:`thewalrus.reference.hafnian()`
 Hafnian of a matrix with repeated rows and columns                               :func:`thewalrus.hafnian_repeated()`
-Hafnians of all possible reductions of a given matrix                             :func:`thewalrus.hafnian_batched()`
+Hafnians of all possible reductions of a given matrix                            :func:`thewalrus.hafnian_batched()`
 Hafnian samples of a Gaussian state                                              :func:`thewalrus.samples.hafnian_sample_state()`
 Torontonian samples of a Gaussian state                                          :func:`thewalrus.samples.torontonian_sample_state()`
 Hafnian samples of a graph                                                       :func:`thewalrus.samples.hafnian_sample_graph()`
@@ -21,6 +21,7 @@ All probability amplitudes of a pure Gaussian state
 All matrix elements of a general Gaussian state                                  :func:`thewalrus.quantum.density_matrix()`
 A particular probability amplitude of pure Gaussian state                        :func:`thewalrus.quantum.pure_state_amplitude()`
 A particular matrix element of a general Gaussian state                          :func:`thewalrus.quantum.density_matrix_element()`
+The Fock representation of a Gaussian unitary                                    :func:`thewalrus.quantum.fock_tensor()`
 ================================================================================ =============================
 
 Note that all the hafnian functions listed above generalize to loop hafnians.

include/hermite_multidimensional.hpp

@@ -168,7 +168,6 @@ inline std::vector<T> hermite_multidimensional_cpp(std::vector<T> &R_mat, std::v
     std::vector<double> factors(resolution+1, 0);
     int jump = 0;
 
-
     for (ullint jj = 0; jj < Hdim-1; jj++) {
 
         if (jump > 0) {
@@ -205,9 +204,7 @@ inline std::vector<T> hermite_multidimensional_cpp(std::vector<T> &R_mat, std::v
         ullint fromCoordinate = vec2index(jumpFrom, resolution);
 
 
-        for (int ii = 0; ii < dim; ii++) {
-            H[nextCoordinate] = H[nextCoordinate] + R(k, ii) * y(ii, 0);
-        }
+	H[nextCoordinate] = H[nextCoordinate] + y(k, 0);
         H[nextCoordinate] = H[nextCoordinate] * H[fromCoordinate];
 
         std::vector<int> tmpjump(dim, 0);

tests/libwalrus_unittest.cpp

@@ -734,7 +734,7 @@ TEST(LoopHafnianEigenComplex, EvenOnes) {
 }
 
 
-// Check loop hafnian with eignevalues for random complex matrices with even dimensions.
+// Check loop hafnian with eigenvalues for random complex matrices with even dimensions.
 TEST(LoopHafnianEigenComplex, EvenRandom) {
     std::vector<std::complex<double>> mat2(4, 0.0);
     std::vector<std::complex<double>> mat4(16, 0.0);
@@ -776,7 +776,7 @@ TEST(LoopHafnianEigenComplex, EvenRandom) {
 
 }
 
-// Check loop hafnian with eignevalues for complex matrices with odd dimensions.
+// Check loop hafnian with eigenvalues for complex matrices with odd dimensions.
 TEST(LoopHafnianEigenComplex, Odd) {
     std::vector<std::complex<double>> mat3(9, std::complex<double>(1.0, 0.0));
     std::vector<std::complex<double>> mat5(25, std::complex<double>(1.0, 0.0));
@@ -1139,7 +1139,7 @@ TEST(TorontonianDouble, Analytical) {
 
 
 namespace batchhafnian {
-TEST(BatchHafnian, CompleteGraph) {
+TEST(BatchHafnian, Clements) {
     std::vector<std::complex<double>> mat4{std::complex<double>(-0.28264629150778969, 0.39867701584672210), std::complex<double>(-0.06086128222348247, -0.12220227033305252), std::complex<double>(-0.22959477315790058, 0.00000000000000008), std::complex<double>(-0.00660678867199307, -0.09884501458235322), std::complex<double>(-0.06086128222348247, -0.12220227033305252), std::complex<double>(0.38245649793510783, -0.41413300040003126), std::complex<double>(-0.00660678867199307, 0.09884501458235322), std::complex<double>(-0.13684045954832844, 0.00000000000000006), std::complex<double>(-0.22959477315790058, -0.00000000000000008), std::complex<double>(-0.00660678867199307, 0.09884501458235322), std::complex<double>(-0.28264629150778969, -0.39867701584672210), std::complex<double>(-0.06086128222348247, 0.12220227033305252), std::complex<double>(-0.00660678867199307, -0.09884501458235322), std::complex<double>(-0.13684045954832844, -0.00000000000000006), std::complex<double>(-0.06086128222348247, +0.12220227033305252), std::complex<double>(0.38245649793510783, 0.41413300040003126)};
     std::vector<std::complex<double>> d4{std::complex<double>(0.66917130190858, -1.52776303400764), std::complex<double>(-2.95847055822102, -1.29582519437023), std::complex<double>(0.66917130190858, 1.52776303400764), std::complex<double>(-2.95847055822102, 1.29582519437023)};
     std::vector<std::complex<double>> out(256, 0.0);
@@ -1273,8 +1273,17 @@ TEST(BatchHafnian, CompleteGraph) {
                                     9.62872951e-01, -2.53217806e+00,  3.10967275e+00, -1.42559575e-14};
     int res = 4;
     int renorm = 0;
-
-    out = libwalrus::hermite_multidimensional_cpp(mat4, d4, res, renorm);
+    std::vector<std::complex<double>> d4c{std::complex<double>(0.0, 0.0), std::complex<double>(0.0, 0.0), std::complex<double>(0.0, 0.0), std::complex<double>(0.0, 0.0)};
+    for(int i=0; i<4; i++) {
+        for(int j=0; j<4; j++) {
+            d4c[i] += mat4[4*i+j] * d4[j];
+        }
+    }
+    // Note that internaly we are implementing the modified multidimensional
+    // Hermite polynomials, which means that we have to mat4 * d4 as the
+    // second argument, this is precisely what is done in the previous two loops
+
+    out = libwalrus::hermite_multidimensional_cpp(mat4, d4c, res, renorm);
 
     for (int i = 0; i < 256; i++) {
         EXPECT_NEAR(expected_re[i], std::real(out[i]), tol2);

thewalrus/_hermite_multidimensional.py

@@ -14,52 +14,13 @@
 """
 Hermite Multidimensional Python interface
 """
-from itertools import product
 import numpy as np
 
 from .libwalrus import hermite_multidimensional as hm
 from ._hafnian import input_validation
 
-
-def return_prod(C, index):
-    r"""Given an array :math:`C_{i,j}` and an array or list of indices
-    :math:`index = [i_1,i_2,i_3,\dots,i_n] `, returns :math:`prod_{k=1}^n C_{k,i_k}`.
-
-    Args:
-        C (array): An array
-        index (array): A set of indices
-
-    Returns:
-        complex: the product of the array elements determined by index
-    """
-    return np.prod([C[mode, val] for mode, val in enumerate(index)])
-
-
-def expansion_coeff(alpha, cutoff, renorm=True):
-    r"""Returns the (quasi) geometric series as a vector with components
-    :math:`\alpha^i/\sqrt{i!}` for :math:`0 \leq i < \texttt{cutoff}`.
-
-    Args:
-        alpha (complex): ratio of the geometric series
-        cutoff (int): cutoff truncation of the geometric series
-        renorm (bool): if ``False``, the components are not normalized
-            by the square-root factorials
-
-    Returns:
-        array: the (quasi) geometric series
-    """
-    vals = np.empty([cutoff], dtype=type(alpha))
-    vals[0] = 1.0
-    if renorm:
-        for i in range(1, cutoff):
-            vals[i] = vals[i - 1] * alpha / np.sqrt(i)
-    else:
-        for i in range(1, cutoff):
-            vals[i] = vals[i - 1] * alpha
-    return vals
-
-
-def hermite_multidimensional(R, cutoff, y=None, renorm=False, make_tensor=True):
+# pylint: disable=too-many-arguments
+def hermite_multidimensional(R, cutoff, y=None, renorm=False, make_tensor=True, modified=False):
     r"""Returns the multidimensional Hermite polynomials :math:`H_k^{(R)}(y)`.
 
     Here :math:`R` is an :math:`n \times n` square matrix, and
@@ -84,14 +45,25 @@ def hermite_multidimensional(R, cutoff, y=None, renorm=False, make_tensor=True):
         y (array): vector argument of the Hermite polynomial
         renorm (bool): If ``True``, normalizes the returned multidimensional Hermite
             polynomials such that :math:`H_k^{(R)}(y)/\prod_i k_i!`
-        make_tensor: If ``False``, returns a flattened one dimensional array
+        make_tensor (bool): If ``False``, returns a flattened one dimensional array
             containing the values of the polynomial
+        modified (bool): whether to return the modified multidimensional Hermite polynomials or the standard ones
 
     Returns:
         (array): the multidimensional Hermite polynomials
     """
+
     input_validation(R)
     n, _ = R.shape
+
+    if (modified is False) and (y is not None):
+        m = y.shape[0]
+        if m == n:
+            ym = R @ y
+            return hermite_multidimensional(
+                R, cutoff, y=ym, renorm=renorm, make_tensor=make_tensor, modified=True
+            )
+
     if y is None:
         y = np.zeros([n], dtype=complex)
 
@@ -107,7 +79,7 @@ def hermite_multidimensional(R, cutoff, y=None, renorm=False, make_tensor=True):
 
     return values
 
-
+# pylint: disable=too-many-arguments
 def hafnian_batched(A, cutoff, mu=None, tol=1e-12, renorm=False, make_tensor=True):
     r"""Calculates the hafnian of :func:`reduction(A, k) <hafnian.reduction>`
     for all possible values of vector ``k`` below the specified cutoff.
@@ -141,36 +113,17 @@ def hafnian_batched(A, cutoff, mu=None, tol=1e-12, renorm=False, make_tensor=Tru
     Returns:
         (array): the values of the hafnians for each value of :math:`k` up to the cutoff
     """
-    # pylint: disable=too-many-return-statements,too-many-branches,too-many-arguments
     input_validation(A, tol=tol)
     n, _ = A.shape
 
-    if not np.allclose(A, np.zeros([n, n])):
-        if mu is not None:
-            try:
-                yi = np.linalg.solve(A, mu)
-            except np.linalg.LinAlgError:
-                raise ValueError("The matrix does not have an inverse")
-            return hermite_multidimensional(
-                -A, cutoff, y=-yi, renorm=renorm, make_tensor=make_tensor
-            )
-        yi = np.zeros([n], dtype=complex)
+    if mu is not None:
         return hermite_multidimensional(
-            -A, cutoff, y=-yi, renorm=renorm, make_tensor=make_tensor
+            -A, cutoff, y=mu, renorm=renorm, make_tensor=make_tensor, modified=True
         )
-    # Note the minus signs in the arguments. Those are intentional and are due to the fact that Dodonov et al. in PRA 50, 813 (1994) use (p,q) ordering instead of (q,p) ordering
-
-    if mu is None:
-        tensor = np.zeros([cutoff ** n], dtype=complex)
-        tensor[0] = 1.0
-    else:
-        index = cutoff * np.ones([n], dtype=int)
-        tensor = np.empty(index, dtype=complex)
-        prim = np.array([expansion_coeff(alpha, cutoff, renorm=renorm) for alpha in mu])
-        for i in product(range(cutoff), repeat=n):
-            tensor[i] = return_prod(prim, i)
+    yi = np.zeros([n], dtype=complex)
+    return hermite_multidimensional(
+        -A, cutoff, y=yi, renorm=renorm, make_tensor=make_tensor, modified=True
+    )
 
-    if make_tensor:
-        return tensor
 
-    return tensor.flatten()
+# Note the minus signs in the arguments. Those are intentional and are due to the fact that Dodonov et al. in PRA 50, 813 (1994) use (p,q) ordering instead of (q,p) ordering

thewalrus/quantum.py

@@ -380,17 +380,12 @@ def density_matrix(mu, cov, post_select=None, normalize=False, cutoff=5, hbar=2)
         sf_order = tuple(chain.from_iterable([[i, i + N] for i in range(N)]))
 
         if np.allclose(mu, np.zeros_like(mu)):
-            tensor = pref * hermite_multidimensional(-A, cutoff, renorm=True)
+            tensor = pref * hermite_multidimensional(-A, cutoff, renorm=True, modified=True)
             return tensor.transpose(sf_order)
-        try:
-            beta = Beta(mu)
-            y = np.linalg.inv(A) @ (beta - A @ beta.conj())
-            tensor = pref * hermite_multidimensional(-A, cutoff, y=-y, renorm=True)
-            return tensor.transpose(sf_order)
-
-        except np.linalg.LinAlgError:
-            pass
-        post_select = {}
+        beta = Beta(mu)
+        y = beta - A @ beta.conj()
+        tensor = pref * hermite_multidimensional(-A, cutoff, y=y, renorm=True, modified=True)
+        return tensor.transpose(sf_order)
 
     M = N - len(post_select)
     rho = np.zeros([cutoff] * (2 * M), dtype=np.complex128)

thewalrus/tests/test_hermite_multidimensional.py

@@ -85,7 +85,6 @@ def test_hafnian_batched_loops():
         ]
     )
     expected = hafnian_batched(A, n_photon, mu=mu, make_tensor=False)
-
     assert np.allclose(expected, v1)
 
 
@@ -121,3 +120,15 @@ def test_hafnian_batched_zero_loops_no_edges():
     expected = hafnian_batched(A, n_photon, make_tensor=False)
 
     assert np.allclose(expected, v1)
+
+
+def test_hermite_vs_hermite_modified():
+    """Test the relation hermite and hermite modified"""
+    n_modes = 2
+    A = np.zeros([n_modes, n_modes], dtype=complex)
+    mu = np.random.rand(n_modes) + 1j * np.random.rand(n_modes)
+    cutoff = 3
+    assert np.allclose(
+        hermite_multidimensional(A, cutoff, y=A @ mu, modified=True),
+        hermite_multidimensional(A, cutoff, y=mu),
+    )