passive transformation

.github/CHANGELOG.md

@@ -2,6 +2,10 @@
 
 ### New features
 
+* ``symplectic.squeezing`` function is now generalized to multiple modes of single mode squeezing [#249](https://github.com/XanaduAI/thewalrus/pull/249)
+* adds a function ``symplectic.passive_transformation`` which allows for Gaussian states to be transformed by arbitrary non-unitary, non-square linear optical transformations [#249](https://github.com/XanaduAI/thewalrus/pull/249)
+* ``torontonian_sample_state`` now can sample displaced Gaussian states [#248](https://github.com/XanaduAI/thewalrus/pull/248)
+
 ### Improvements
 
 * Speeds up the calculation of photon number variances/covariances [#224](https://github.com/XanaduAI/thewalrus/pull/224)
@@ -14,7 +18,7 @@
 
 This release contains contributions from (in alphabetical order):
 
-Nicolas Quesada
+Jake Bulmer, Nicolas Quesada
 
 ---
 

thewalrus/symplectic.py

@@ -45,6 +45,7 @@
 
 .. autosummary::
     two_mode_squeezing
+    squeezing
     interferometer
     loss
     mean_photon_number
@@ -99,8 +100,33 @@ def expand_vector(alpha, mode, N, hbar=2.0):
     return r
 
 
+def expand_passive(T, modes, N):
+    r"""Returns and expanded linear optical transformation
+    acting on specified modes, with identity acting on all other modes
+
+    Args:
+        T (array): square :math:`M \times M` matrix of linear optical transformation
+        modes (array): :math:`M` modes of the tranformation
+        N (int): number of modes
+
+    Returns:
+        array: phase-space displacement vector of size 2*N
+    """
+
+    if T.shape[0] != T.shape[1]:
+        raise ValueError("The input matrix is not square")
+
+    if len(modes) != T.shape[0]:
+        raise ValueError("length of modes must match the shape of T")
+
+    T_expand = np.eye(N, dtype=T.dtype)
+    T_expand[np.ix_(modes, modes)] = T
+
+    return T_expand
+
+
 def reduced_state(mu, cov, modes):
-    r""" Returns the vector of means and the covariance matrix of the specified modes.
+    r"""Returns the vector of means and the covariance matrix of the specified modes.
 
     Args:
         modes (int of Sequence[int]): indices of the requested modes
@@ -146,22 +172,35 @@ def vacuum_state(modes, hbar=2.0, dtype=np.float64):
     return state
 
 
-def squeezing(r, phi, dtype=np.float64):
+def squeezing(r, phi=None, dtype=np.float64):
     r"""Squeezing. In fock space this corresponds to \exp(\tfrac{1}{2}r e^{i \phi} (a^2 - a^{\dagger 2}) ).
 
     Args:
-        r (float): squeezing magnitude
-        phi (float): rotation parameter
+        r (array): squeezing magnitude
+        phi (array): rotation parameter
         dtype (numpy.dtype): datatype to represent the Symplectic matrix
     Returns:
         array: symplectic transformation matrix
     """
     # pylint: disable=assignment-from-no-return
-    cp = np.cos(phi, dtype=dtype)
-    sp = np.sin(phi, dtype=dtype)
-    ch = np.cosh(r, dtype=dtype)
-    sh = np.sinh(r, dtype=dtype)
-    S = np.array([[ch - cp * sh, -sp * sh], [-sp * sh, ch + cp * sh]])
+
+    r = np.atleast_1d(r)
+
+    if phi is None:
+        phi = np.zeros_like(r)
+
+    phi = np.atleast_1d(phi)
+
+    M = len(r)
+    S = np.identity(2 * M, dtype=dtype)
+
+    for i in range(M):
+        r_i = r[i]
+        phi_i = phi[i]
+        S[i, i] = np.cosh(r_i) - np.sinh(r_i) * np.cos(phi_i)
+        S[i, i + M] = -np.sinh(r_i) * np.sin(phi_i)
+        S[i + M, i] = -np.sinh(r_i) * np.sin(phi_i)
+        S[i + M, i + M] = np.cosh(r_i) + np.sinh(r_i) * np.cos(phi_i)
 
     return S
 
@@ -211,6 +250,28 @@ def interferometer(U):
     return S
 
 
+def passive_transformation(mu, cov, T, hbar=2):
+    r"""
+    covariance matrix transformation for arbitrary linear optical channel
+    mapping an :math:`N` modes state to an :math:`M` mode state
+    Args:
+        mu (array): :math:`2N`-length means vector
+        cov (array): :math:`2N \times 2N` covariance matrix
+        T (array): :math:`M \times N` linear optical transformation
+    Returns:
+        mu (array): :math:`2M`-length transformed means vector
+        cov (array): :math:`2M \times 2M` tranformed covariance matrix
+    """
+
+    P = interferometer(T)
+    L = (hbar / 2) * (np.eye(P.shape[0]) - P @ P.T)
+
+    cov = P @ cov @ P.T + L
+    mu = P @ mu
+
+    return mu, cov
+
+
 # pylint: disable=too-many-arguments
 def loss(mu, cov, T, mode, nbar=0, hbar=2):
     r"""Loss channel acting on a Gaussian state.
@@ -241,6 +302,7 @@ def loss(mu, cov, T, mode, nbar=0, hbar=2):
 
     return mu_res, cov_res
 
+
 ### Comment: This function strongly overlaps with `quantum.photon_number_mean`
 ### Wonder if it is worth removing it.
 def mean_photon_number(mu, cov, hbar=2):
@@ -314,7 +376,7 @@ def sympmat(N, dtype=np.float64):
 
 
 def is_symplectic(S, rtol=1e-05, atol=1e-08):
-    r""" Checks if matrix S is a symplectic matrix
+    r"""Checks if matrix S is a symplectic matrix
 
     Args:
         S (array): a square matrix
@@ -332,6 +394,7 @@ def is_symplectic(S, rtol=1e-05, atol=1e-08):
 
     return np.allclose(S.T @ Omega @ S, Omega, rtol=rtol, atol=atol)
 
+
 def autonne(A, rtol=1e-05, atol=1e-08, svd_order=True):
     r"""Autonne-Takagi decomposition of a complex symmetric (not Hermitian!) matrix.
 
@@ -364,6 +427,7 @@ def autonne(A, rtol=1e-05, atol=1e-08, svd_order=True):
         return (vals[n : 2 * n])[::-1], U[:, ::-1]
     return vals[n : 2 * n], U
 
+
 def xxpp_to_xpxp(S):
     """Permutes the entries of the input from xxpp ordering to xpxp ordering.
 
@@ -389,6 +453,7 @@ def xxpp_to_xpxp(S):
 
     return S[ind]
 
+
 def xpxp_to_xxpp(S):
     """Permutes the entries of the input from xpxp ordering to xxpp ordering.
 

thewalrus/tests/test_symplectic.py

@@ -19,6 +19,7 @@
 from scipy.linalg import block_diag
 
 from thewalrus import symplectic
+from thewalrus.quantum import is_valid_cov
 
 
 # pylint: disable=too-few-public-methods
@@ -59,6 +60,28 @@ def test_squeezing(self, tol):
         expected = rotation @ np.diag(np.exp([-2*r, 2*r])) @ rotation.T
         assert np.allclose(out, expected, atol=tol, rtol=0)
 
+    def test_squeezing_no_phi(self, tol):
+        """Test the squeezing symplectic transform without specifying phi"""
+        r = 0.543
+        phi = 0.
+        S = symplectic.squeezing(r)
+        out = S @ S.T
+
+        # apply to an identity covariance matrix
+        rotation = np.array([[np.cos(phi/2), -np.sin(phi/2)], [np.sin(phi/2), np.cos(phi/2)]])
+        expected = rotation @ np.diag(np.exp([-2*r, 2*r])) @ rotation.T
+        assert np.allclose(out, expected, atol=tol, rtol=0)
+
+    def test_symplectic_multimode(self, tol):
+        """Test multimode version gives symplectic matrix"""
+        r = [0.543] * 4
+        phi = [0.123] * 4
+        S = symplectic.squeezing(r, phi)
+
+        # the symplectic matrix
+        O = symplectic.sympmat(4)
+
+        assert np.allclose(S @ O @ S.T, O, atol=tol, rtol=0)
 
 class TestTwoModeSqueezing:
     """Tests for the TMS symplectic"""
@@ -176,6 +199,76 @@ def test_rotation(self, tol):
         np.allclose(S, expected, atol=tol, rtol=0)
 
 
+class TestPassiveTransformation:
+    """ tests for linear transformation """
+    def test_transformation(self, tol):
+        """Test that an transformation returns the correct state"""
+
+        M = 4
+        cov = np.arange(4 * M ** 2, dtype=np.float64).reshape((2*M, 2*M))
+        mu = np.arange(2 * M, dtype=np.float64)
+
+        T = np.sqrt(0.9) * M ** (-0.5) * np.ones((6,M), dtype=np.float64)
+
+        mu_out, cov_out = symplectic.passive_transformation(mu, cov, T)
+        # fmt:off
+        expected_mu = np.array([ 2.84604989,  2.84604989,  2.84604989,  2.84604989,  2.84604989,
+                                 2.84604989, 10.43551628, 10.43551628, 10.43551628, 10.43551628,
+                                 10.43551628, 10.43551628])
+        expected_cov = np.array([
+            [ 48.7,  47.7,  47.7,  47.7,  47.7,  47.7,  63. ,  63. ,  63. , 63. ,  63. ,  63. ],
+            [ 47.7,  48.7,  47.7,  47.7,  47.7,  47.7,  63. ,  63. ,  63. , 63. ,  63. ,  63. ],
+            [ 47.7,  47.7,  48.7,  47.7,  47.7,  47.7,  63. ,  63. ,  63. , 63. ,  63. ,  63. ],
+            [ 47.7,  47.7,  47.7,  48.7,  47.7,  47.7,  63. ,  63. ,  63. , 63. ,  63. ,  63. ],
+            [ 47.7,  47.7,  47.7,  47.7,  48.7,  47.7,  63. ,  63. ,  63. , 63. ,  63. ,  63. ],
+            [ 47.7,  47.7,  47.7,  47.7,  47.7,  48.7,  63. ,  63. ,  63. , 63. ,  63. ,  63. ],
+            [163.8, 163.8, 163.8, 163.8, 163.8, 163.8, 178.3, 177.3, 177.3, 177.3, 177.3, 177.3],
+            [163.8, 163.8, 163.8, 163.8, 163.8, 163.8, 177.3, 178.3, 177.3, 177.3, 177.3, 177.3],
+            [163.8, 163.8, 163.8, 163.8, 163.8, 163.8, 177.3, 177.3, 178.3, 177.3, 177.3, 177.3],
+            [163.8, 163.8, 163.8, 163.8, 163.8, 163.8, 177.3, 177.3, 177.3, 178.3, 177.3, 177.3],
+            [163.8, 163.8, 163.8, 163.8, 163.8, 163.8, 177.3, 177.3, 177.3, 177.3, 178.3, 177.3],
+            [163.8, 163.8, 163.8, 163.8, 163.8, 163.8, 177.3, 177.3, 177.3, 177.3, 177.3, 178.3]])
+        # fmt:on
+
+        assert np.allclose(mu_out, expected_mu, atol=tol, rtol=0)
+        assert np.allclose(cov_out, expected_cov, atol=tol, rtol=0)
+
+    @pytest.mark.parametrize("M", range(1,10))
+    def test_valid_cov(self, M, tol):
+        """test that the output is a valid covariance matrix, even when not square"""
+        a = np.arange(4 * M ** 2, dtype=np.float64).reshape((2*M, 2*M))
+        cov = a @ a.T + np.eye(2*M)
+        mu = np.arange(2 * M, dtype=np.float64)
+
+        T = np.sqrt(0.9) * M ** (-0.5) * np.ones((6,M), dtype=np.float64)
+
+        mu_out, cov_out = symplectic.passive_transformation(mu, cov, T)
+
+        assert cov_out.shape == (12, 12)
+        assert len(mu_out) == 12
+        assert is_valid_cov(cov_out, atol=tol, rtol=0)
+
+    @pytest.mark.parametrize("M", range(1, 6))
+    def test_unitary(self, M, tol):
+        """
+        test that the outputs agree with the interferometer class when
+        transformation is unitary
+        """
+        a = np.arange(4 * M ** 2, dtype=np.float64).reshape((2 * M, 2 * M))
+        cov = a @ a.T + np.eye(2 * M)
+        mu = np.arange(2 * M, dtype=np.float64)
+
+        U = M ** (-0.5) * np.fft.fft(np.eye(M))
+        S_U = symplectic.interferometer(U)
+        cov_U = S_U @ cov @ S_U.T
+        mu_U = S_U @ mu
+
+        mu_T, cov_T = symplectic.passive_transformation(mu, cov, U)
+
+        assert np.allclose(mu_U, mu_T, atol=tol, rtol=0)
+        assert np.allclose(cov_U, cov_T, atol=tol, rtol=0)
+
+
 class TestReducedState:
     """Tests for the reduced state function"""
 
@@ -465,6 +558,30 @@ def test_expand_one(self, hbar, tol):
         expected[mode + N] = np.sqrt(2 * hbar) * alpha.imag
         assert np.allclose(r, expected, atol=tol, rtol=0)
 
+class TestExpandPassive:
+    """Tests for expanding a displacement operation into a phase-space displacement vector"""
+
+    def test_expand_one(self, tol):
+        """Test that a 1x1 matrix is expanded correctly"""
+        T = np.array([[0.5]])
+
+        T_expand = symplectic.expand_passive(T, [1], 3)
+
+        expected = np.array([[1,0,0],[0,0.5,0],[0,0,1]])
+
+        assert np.allclose(T_expand, expected, atol=tol, rtol=0)
+
+    def test_expend_not_square(self):
+        """ test that error is raised for non square input"""
+        with pytest.raises(ValueError, match="The input matrix is not square"):
+            symplectic.expand_passive(np.ones((3,2)), [0,1,2], 5)
+
+    def test_modes_length(self):
+        """ test that error is raised when length of modes array is incorrect"""
+        with pytest.raises(ValueError, match="length of modes must match the shape of T"):
+            symplectic.expand_passive(np.ones((3,3)), [0,1,2,3,4], 8)
+
+
 class TestSymplecticExpansion:
     """Tests for the expanding a symplectic matrix"""