diff --git a/docs/quick_guide.rst b/docs/quick_guide.rst index 9b3272593..67d297c19 100644 --- a/docs/quick_guide.rst +++ b/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. diff --git a/include/hermite_multidimensional.hpp b/include/hermite_multidimensional.hpp index a6a560f37..5f95925f6 100644 --- a/include/hermite_multidimensional.hpp +++ b/include/hermite_multidimensional.hpp @@ -168,7 +168,6 @@ inline std::vector hermite_multidimensional_cpp(std::vector &R_mat, std::v std::vector 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 hermite_multidimensional_cpp(std::vector &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 tmpjump(dim, 0); diff --git a/tests/libwalrus_unittest.cpp b/tests/libwalrus_unittest.cpp index c84584884..4c24a6697 100644 --- a/tests/libwalrus_unittest.cpp +++ b/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> mat2(4, 0.0); std::vector> 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> mat3(9, std::complex(1.0, 0.0)); std::vector> mat5(25, std::complex(1.0, 0.0)); @@ -1139,7 +1139,7 @@ TEST(TorontonianDouble, Analytical) { namespace batchhafnian { -TEST(BatchHafnian, CompleteGraph) { +TEST(BatchHafnian, Clements) { std::vector> mat4{std::complex(-0.28264629150778969, 0.39867701584672210), std::complex(-0.06086128222348247, -0.12220227033305252), std::complex(-0.22959477315790058, 0.00000000000000008), std::complex(-0.00660678867199307, -0.09884501458235322), std::complex(-0.06086128222348247, -0.12220227033305252), std::complex(0.38245649793510783, -0.41413300040003126), std::complex(-0.00660678867199307, 0.09884501458235322), std::complex(-0.13684045954832844, 0.00000000000000006), std::complex(-0.22959477315790058, -0.00000000000000008), std::complex(-0.00660678867199307, 0.09884501458235322), std::complex(-0.28264629150778969, -0.39867701584672210), std::complex(-0.06086128222348247, 0.12220227033305252), std::complex(-0.00660678867199307, -0.09884501458235322), std::complex(-0.13684045954832844, -0.00000000000000006), std::complex(-0.06086128222348247, +0.12220227033305252), std::complex(0.38245649793510783, 0.41413300040003126)}; std::vector> d4{std::complex(0.66917130190858, -1.52776303400764), std::complex(-2.95847055822102, -1.29582519437023), std::complex(0.66917130190858, 1.52776303400764), std::complex(-2.95847055822102, 1.29582519437023)}; std::vector> 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> d4c{std::complex(0.0, 0.0), std::complex(0.0, 0.0), std::complex(0.0, 0.0), std::complex(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); diff --git a/thewalrus/_hermite_multidimensional.py b/thewalrus/_hermite_multidimensional.py index def31dfab..5faf8d8da 100644 --- a/thewalrus/_hermite_multidimensional.py +++ b/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) ` 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 diff --git a/thewalrus/quantum.py b/thewalrus/quantum.py index 027aae494..28cfb575b 100644 --- a/thewalrus/quantum.py +++ b/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) diff --git a/thewalrus/tests/test_hermite_multidimensional.py b/thewalrus/tests/test_hermite_multidimensional.py index f77ec5e8d..79b1cf329 100644 --- a/thewalrus/tests/test_hermite_multidimensional.py +++ b/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), + )