Source code for thewalrus.symplectic

# Copyright 2019 Xanadu Quantum Technologies Inc.

# 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.
"""
Symplectic operations
=====================

**Module name:** :mod:`thewalrus.symplectic`

.. currentmodule:: thewalrus.symplectic

Contains some Gaussian operations and auxiliary functions.

Auxiliary functions
-------------------

.. autosummary::
    expand
    expand_vector
    expand_passive
    reduced_state
    is_symplectic
    sympmat
    xxpp_to_xpxp
    xpxp_to_xxpp

Gaussian states
---------------

.. autosummary::
    vacuum_state


Gates and operations
--------------------

.. autosummary::
    two_mode_squeezing
    squeezing
    interferometer
    loss
    mean_photon_number
    beam_splitter
    rotation

Code details
------------
"""
import warnings
import numpy as np
from scipy.sparse import (
    identity as sparse_identity,
    issparse,
    coo_array,
    dia_array,
    bsr_array,
    csr_array,
)


[docs] def expand(S, modes, N): r"""Expands a Symplectic matrix S to act on the entire subsystem. If the input is a single mode symplectic, then extends it to act on multiple modes. Supports scipy sparse matrices. Instances of ``coo_array``, ``dia_array``, ``bsr_array`` will be transformed into `csr_array``. Args: S (ndarray or spmatrix): a :math:`2M\times 2M` Symplectic matrix modes (Sequence[int]): the list of modes S acts on N (int): full size of the subsystem Returns: array: the resulting :math:`2N\times 2N` Symplectic matrix """ M = S.shape[0] // 2 S2 = ( np.identity(2 * N, dtype=S.dtype) if not issparse(S) else sparse_identity(2 * N, dtype=S.dtype, format="csr") ) if issparse(S) and isinstance(S, (coo_array, dia_array, bsr_array)): # cast to sparse matrix that supports slicing and indexing warnings.warn( "Unsupported sparse matrix type, returning a Compressed Sparse Row (CSR) matrix." ) S = csr_array(S) w = np.array([modes]) if isinstance(modes, int) else np.array(modes) # extend single mode symplectic to act on selected modes if M == 1: for m in w: S2[m, m], S2[m + N, m + N] = S[0, 0], S[1, 1] # X, P S2[m, m + N], S2[m + N, m] = S[0, 1], S[1, 0] # XP, PX return S2 # make symplectic act on the selected subsystems S2[w.reshape(-1, 1), w.reshape(1, -1)] = S[:M, :M].copy() # X S2[(w + N).reshape(-1, 1), (w + N).reshape(1, -1)] = S[M:, M:].copy() # P S2[w.reshape(-1, 1), (w + N).reshape(1, -1)] = S[:M, M:].copy() # XP S2[(w + N).reshape(-1, 1), w.reshape(1, -1)] = S[M:, :M].copy() # PX return S2
[docs] def expand_vector(alpha, mode, N, hbar=2.0): """Returns the phase-space displacement vector associated to a displacement. Args: alpha (complex): complex displacement mode (int): mode index N (int): number of modes Returns: array: phase-space displacement vector of size 2*N """ alpharealdtype = np.dtype(type(alpha)) r = np.zeros(2 * N, dtype=alpharealdtype) r[mode] = np.sqrt(2 * hbar) * alpha.real r[N + mode] = np.sqrt(2 * hbar) * alpha.imag return r
[docs] def expand_passive(T, modes, N): r"""Returns the 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): the :math:`M` modes of the transformation N (int): number of modes in the new expanded transformation Returns: array: :math:`N \times N` array of expanded passive transformation """ 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
[docs] def reduced_state(mu, cov, 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 Returns: tuple (means, cov): means is an array containing the vector of means, and cov is a square array containing the covariance matrix. Both use the :math:`xp`-ordering. """ N = len(mu) // 2 if modes == list(range(N)): # reduced state is full state return mu, cov # reduce state down to specified subsystems if isinstance(modes, int): modes = [modes] if np.any(np.array(modes) > N): raise ValueError("The specified modes cannot be larger than number of subsystems.") ind = np.concatenate([np.array(modes), np.array(modes) + N]) rows = ind.reshape(-1, 1) cols = ind.reshape(1, -1) return mu.copy()[ind], cov.copy()[rows, cols]
[docs] def vacuum_state(modes, hbar=2.0, dtype=np.float64): r"""Returns the vacuum state. Args: modes (str): Returns the vector of means and the covariance matrix hbar (float): (default 2) the value of :math:`\hbar` in the commutation relation :math:`[\x,\p]=i\hbar` dtype (numpy.dtype): datatype to represent the covariance matrix and vector of means Returns: list[array]: the means vector and covariance matrix of the vacuum state """ means = np.zeros((2 * modes), dtype=dtype) cov = np.identity(2 * modes, dtype=dtype) * hbar / 2 state = [means, cov] return state
[docs] def squeezing(r, phi=None, dtype=np.float64): r"""Squeezing. In fock space this corresponds to: .. math:: \exp(\tfrac{1}{2}r e^{i \phi} (a^2 - a^{\dagger 2}) ). By passing an array of squeezing parameters and phases, it applies a tensor product of squeezing operations. Args: r (Union[array, float]): squeezing magnitude phi (Union[array, float]): rotation parameter. If ``None``, then the function uses zeros of the same shape as ``r``. dtype (numpy.dtype): datatype to represent the Symplectic matrix. Defaults to ``numpy.float64``. Returns: array: symplectic transformation matrix """ # pylint: disable=assignment-from-no-return 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, (r_i, phi_i) in enumerate(zip(r, phi)): 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
[docs] def two_mode_squeezing(r, phi, dtype=np.float64): """Two-mode squeezing. Args: r (float): squeezing magnitude phi (float): 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, 0, sp * sh], [cp * sh, ch, sp * sh, 0], [0, sp * sh, ch, -cp * sh], [sp * sh, 0, -cp * sh, ch], ] ) return S
[docs] def interferometer(U): """Interferometer. Args: U (array): unitary matrix Returns: array: symplectic transformation matrix """ X = U.real Y = U.imag S = np.block([[X, -Y], [Y, X]]) return S
[docs] def passive_transformation(mu, cov, T, hbar=2): r"""Perform a covariance matrix transformation for an arbitrary linear optical channel on an :math:`N` modes state mapping it to a to an :math:`M` modes 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 Keyword Args: hbar (float)=2: the value to use for hbar Returns: array: :math:`2M`-length transformed means vector 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
[docs] def loss(mu, cov, T, mode, nbar=0, hbar=2): r"""Loss channel acting on a Gaussian state. Args: mu (array): means vector cov (array): covariance matri T (float): transmission; 1 corresponds to no loss, 0 to full loss. mode (int): mode to act on nbar (float): thermal mean population (default 0) hbar (float): (default 2) the value of :math:`\hbar` in the commutation relation :math:`[\x,\p]=i\hbar` Returns: tuple[array]: the means vector and covariance matrix of the resulting state """ N = len(cov) // 2 cov_res = cov.copy() mu_res = mu.copy() for m in (mode, mode + N): mu_res[m] *= np.sqrt(T) cov_res[m, :] *= np.sqrt(T) cov_res[:, m] *= np.sqrt(T) cov_res[m, m] += (1 - T) * (2 * nbar + 1) * hbar / 2 return mu_res, cov_res
### Comment: This function strongly overlaps with `quantum.photon_number_mean` ### Wonder if it is worth removing it.
[docs] def mean_photon_number(mu, cov, hbar=2): r"""Calculates the mean photon number for a given one-mode state. Args: mu (array): length-2 vector of means cov (array): :math:`2\times 2` covariance matrix hbar (float): (default 2) the value of :math:`\hbar` in the commutation relation :math:`[\x,\p]=i\hbar` Returns: tuple: the photon number expectation and variance """ ex = (np.trace(cov) + mu.T @ mu) / (2 * hbar) - 1 / 2 var = (np.trace(cov @ cov) + 2 * mu.T @ cov @ mu) / (2 * hbar**2) - 1 / 4 return ex, var
[docs] def beam_splitter(theta, phi, dtype=np.float64): """Beam-splitter. Args: theta (float): transmissivity parameter phi (float): phase parameter dtype (numpy.dtype): datatype to represent the Symplectic matrix Returns: array: symplectic-orthogonal transformation matrix of an interferometer with angles theta and phi """ ct = np.cos(theta, dtype=dtype) st = np.sin(theta, dtype=dtype) eip = np.cos(phi, dtype=dtype) + 1j * np.sin(phi, dtype=dtype) U = np.array( [ [ct, -eip.conj() * st], [eip * st, ct], ] ) return interferometer(U)
[docs] def rotation(theta, dtype=np.float64): """Rotation gate. Args: theta (float): rotation angle dtype (numpy.dtype): datatype to represent the Symplectic matrix Returns: array: rotation matrix by angle theta """ V = np.identity(1) * (np.cos(theta, dtype=dtype) + 1j * np.sin(theta, dtype=dtype)) return interferometer(V)
[docs] def sympmat(N, dtype=np.float64): r"""Returns the matrix :math:`\Omega_n = \begin{bmatrix}0 & I_n\\ -I_n & 0\end{bmatrix}` Args: N (int): positive integer dtype (numpy.dtype): datatype to represent the Symplectic matrix Returns: array: :math:`2N\times 2N` array """ I = np.identity(N, dtype=dtype) O = np.zeros_like(I, dtype=dtype) S = np.block([[O, I], [-I, O]]) return S
[docs] def is_symplectic(S, rtol=1e-05, atol=1e-08): r"""Checks if matrix S is a symplectic matrix Args: S (array): a square matrix Returns: (bool): whether the given matrix is symplectic """ n, m = S.shape if n != m: return False if n % 2 != 0: return False nmodes = n // 2 Omega = sympmat(nmodes, dtype=S.dtype) return np.allclose(S.T @ Omega @ S, Omega, rtol=rtol, atol=atol)
[docs] def xxpp_to_xpxp(S): """Permutes the entries of the input from xxpp ordering to xpxp ordering. Args: S (array): input even dimensional square matrix or array Returns: (array): permuted matrix or array """ shape = S.shape n = shape[0] if n % 2 != 0: raise ValueError("The input array is not even-dimensional") n = n // 2 ind = np.arange(2 * n).reshape(2, -1).T.flatten() if len(shape) == 2: if shape[0] != shape[1]: raise ValueError("The input matrix is not square") return S[:, ind][ind] return S[ind]
[docs] def xpxp_to_xxpp(S): """Permutes the entries of the input from xpxp ordering to xxpp ordering. Args: S (array): input even dimensional square matrix or vector Returns: (array): permuted matrix or vector """ shape = S.shape n = shape[0] if n % 2 != 0: raise ValueError("The input array is not even-dimensional") n = n // 2 ind = np.arange(2 * n).reshape(-1, 2).T.flatten() if len(shape) == 2: if shape[0] != shape[1]: raise ValueError("The input matrix is not square") return S[:, ind][ind] return S[ind]

Contents

Getting started

Background

The Walrus API