Source code for pennylane.devices.default_gaussian
# Copyright 2018-2021 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.# pylint: disable=inconsistent-return-statements"""The :code:`default.gaussian` device is a simulator for Gaussian continuous-variablequantum computations, and can be used as a template for writing PennyLanedevices for new CV backends.It implements the necessary :class:`~pennylane.devices._legacy_device.Device` methods as well as all built-in:mod:`continuous-variable Gaussian operations <pennylane.ops.cv>`, and provides a very simple simulation of aGaussian-based quantum circuit architecture."""importcmath# pylint: disable=attribute-defined-outside-init,too-many-argumentsimportmathimportnumpyasnpfromscipy.specialimportfactorialasfacimportpennylaneasqmlfrompennylane.opsimportIdentityfrom.._versionimport__version__from._legacy_deviceimportDevice# tolerance for numerical errorstolerance=1e-10# ========================================================# auxillary functions# ========================================================
[docs]defpartitions(s,include_singles=True):"""Partitions a sequence into all groupings of pairs and singles of elements. Args: s (sequence): the sequence to partition include_singles (bool): if False, only partitions into pairs is returned. Returns: tuple: returns a nested tuple, containing all partitions of the sequence. """# pylint: disable=too-many-branchesiflen(s)==2:ifinclude_singles:yield(s[0],),(s[1],)yield(tuple(s),)else:# pull off a single item and partition the restifinclude_singles:iflen(s)>1:item_partition=(s[0],)rest=s[1:]rest_partitions=partitions(rest,include_singles)forpinrest_partitions:yield((item_partition),)+pelse:yield(tuple(s),)# pull off a pair of items and partition the restforidx1inrange(1,len(s)):item_partition=(s[0],s[idx1])rest=s[1:idx1]+s[idx1+1:]rest_partitions=partitions(rest,include_singles)forpinrest_partitions:yield((item_partition),)+p
[docs]deffock_prob(cov,mu,event,hbar=2.0):r"""Returns the probability of detection of a particular PNR detection event. For more details, see: * Kruse, R., Hamilton, C. S., Sansoni, L., Barkhofen, S., Silberhorn, C., & Jex, I. "A detailed study of Gaussian Boson Sampling." `arXiv:1801.07488. (2018). <https://arxiv.org/abs/1801.07488>`_ * Hamilton, C. S., Kruse, R., Sansoni, L., Barkhofen, S., Silberhorn, C., & Jex, I. "Gaussian boson sampling." `Physical review letters, 119(17), 170501. (2017). <https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.119.170501>`_ Args: cov (array): :math:`2N\times 2N` covariance matrix mu (array): length-:math:`2N` means vector event (array): length-:math:`N` array of non-negative integers representing the PNR detection event of the multi-mode system. hbar (float): (default 2) the value of :math:`\hbar` in the commutation relation :math:`[\x,\p]=i\hbar`. Returns: float: probability of detecting the event """# number of modesN=len(mu)//2I=np.identity(N)# mean displacement of each modealpha=(mu[:N]+1j*mu[N:])/math.sqrt(2*hbar)# the expectation values (<a_1>, <a_2>,...,<a_N>, <a^\dagger_1>, ..., <a^\dagger_N>)beta=np.concatenate([alpha,alpha.conj()])x=cov[:N,:N]*2/hbarxp=cov[:N,N:]*2/hbarp=cov[N:,N:]*2/hbar# the (Hermitian) matrix elements <a_i^\dagger a_j>aidaj=(x+p+1j*(xp-xp.T)-2*I)/4# the (symmetric) matrix elements <a_i a_j>aiaj=(x-p+1j*(xp+xp.T))/4# calculate the covariance matrix sigma_Q appearing in the Q function:# Q(alpha) = exp[-(alpha-beta).sigma_Q^{-1}.(alpha-beta)/2]/|sigma_Q|Q=np.block([[aidaj,aiaj.conj()],[aiaj,aidaj.conj()]])+np.identity(2*N)# inverse Q matrixQinv=np.linalg.inv(Q)# 1/sqrt(|Q|)sqrt_Qdet=1/math.sqrt(np.linalg.det(Q).real)prefactor=cmath.exp(-beta@Qinv@beta.conj()/2)ifnp.all(np.array(event)==0):# all PNRs detect the vacuum statereturn(prefactor*sqrt_Qdet).real/np.prod(fac(event))# the matrix X_n = [[0, I_n], [I_n, 0]]O=np.zeros_like(I)X=np.block([[O,I],[I,O]])gamma=X@Qinv.conj()@beta# For each mode, repeat the mode number event[i] timesind=[iforsublistin[[idx]*jforidx,jinenumerate(event)]foriinsublist]# extend the indices for xp-ordering of the Gaussian stateind+=[i+Nforiinind]ifnp.linalg.norm(beta)<tolerance:# state has no displacementpart=partitions(ind,include_singles=False)else:part=partitions(ind,include_singles=True)# calculate Hamilton's A matrix: A = X.(I-Q^{-1})*A=X@(np.identity(2*N)-Qinv).conj()summation=np.sum([np.prod([gamma[i[0]]iflen(i)==1elseA[i]foriinp])forpinpart])return(prefactor*sqrt_Qdet*summation).real/np.prod(fac(event))
[docs]defrotation(phi):"""Rotation in the phase space. Args: phi (float): rotation parameter Returns: array: symplectic transformation matrix """returnnp.array([[math.cos(phi),-math.sin(phi)],[math.sin(phi),math.cos(phi)]])
[docs]defdisplacement(state,wire,alpha,hbar=2):"""Displacement in the phase space. Args: state (tuple): contains covariance matrix and means vector wire (int): wire that the displacement acts on alpha (float): complex displacement Returns: tuple: contains the covariance matrix and the vector of means """mu=state[1]mu[wire]+=alpha.real*math.sqrt(2*hbar)mu[wire+len(mu)//2]+=alpha.imag*math.sqrt(2*hbar)returnstate[0],mu
[docs]defsqueezing(r,phi):"""Squeezing in the phase space. Args: r (float): squeezing magnitude phi (float): rotation parameter Returns: array: symplectic transformation matrix """cp=math.cos(phi)sp=math.sin(phi)ch=math.cosh(r)sh=math.sinh(r)returnnp.array([[ch-cp*sh,-sp*sh],[-sp*sh,ch+cp*sh]])
# ========================================================# Arbitrary states and operators# ========================================================
[docs]defsqueezed_cov(r,phi,hbar=2):r"""Returns the squeezed covariance matrix of a squeezed state. Args: r (float): the squeezing magnitude p (float): the squeezing phase :math:`\phi` hbar (float): (default 2) the value of :math:`\hbar` in the commutation relation :math:`[\x,\p]=i\hbar` Returns: array: the squeezed state """cov=np.array([[math.exp(-2*r),0],[0,math.exp(2*r)]])*hbar/2R=rotation(phi/2)returnR@cov@R.T
[docs]defvacuum_state(wires,hbar=2.0):r"""Returns the vacuum state. Args: wires (int): the number of wires to initialize in the vacuum state hbar (float): (default 2) the value of :math:`\hbar` in the commutation relation :math:`[\x,\p]=i\hbar` Returns: array: the vacuum state """means=np.zeros((2*wires))cov=np.identity(2*wires)*hbar/2state=[cov,means]returnstate
[docs]defcoherent_state(a,phi=0,hbar=2.0):r"""Returns a coherent state. Args: a (complex) : the displacement phi (float): the phase hbar (float): (default 2) the value of :math:`\hbar` in the commutation relation :math:`[\x,\p]=i\hbar` Returns: array: the coherent state """alpha=a*cmath.exp(1j*phi)means=np.array([alpha.real,alpha.imag])*math.sqrt(2*hbar)cov=np.identity(2)*hbar/2state=[cov,means]returnstate
[docs]defsqueezed_state(r,phi,hbar=2.0):r"""Returns a squeezed state. Args: r (float): the squeezing magnitude phi (float): the squeezing phase :math:`\phi` hbar (float): (default 2) the value of :math:`\hbar` in the commutation relation :math:`[\x,\p]=i\hbar` Returns: array: the squeezed state """means=np.zeros((2))state=[squeezed_cov(r,phi,hbar),means]returnstate
[docs]defdisplaced_squeezed_state(a,phi_a,r,phi_r,hbar=2.0):r"""Returns a squeezed coherent state Args: a (real): the displacement magnitude phi_a (real): the displacement phase r (float): the squeezing magnitude phi_r (float): the squeezing phase :math:`\phi_r` hbar (float): (default 2) the value of :math:`\hbar` in the commutation relation :math:`[\x,\p]=i\hbar` Returns: array: the squeezed coherent state """alpha=a*cmath.exp(1j*phi_a)means=np.array([alpha.real,alpha.imag])*math.sqrt(2*hbar)state=[squeezed_cov(r,phi_r,hbar),means]returnstate
[docs]defthermal_state(nbar,hbar=2.0):r"""Returns a thermal state. Args: nbar (float): the mean photon number hbar (float): (default 2) the value of :math:`\hbar` in the commutation relation :math:`[\x,\p]=i\hbar` Returns: array: the thermal state """means=np.zeros([2])state=[(2*nbar+1)*np.identity(2)*hbar/2,means]returnstate
[docs]defgaussian_state(cov,mu,hbar=2.0):r"""Returns a Gaussian state. This is simply a bare wrapper function, since the covariance matrix and means vector can be passed via the parameters unchanged. Note that both the covariance and means vector matrix should be in :math:`(\x_1,\dots, \x_N, \p_1, \dots, \p_N)` ordering. Args: cov (array): covariance matrix. Must be dimension :math:`2N\times 2N`, where N is the number of modes mu (array): vector means. Must be length-:math:`2N`, where N is the number of modes hbar (float): (default 2) the value of :math:`\hbar` in the commutation relation :math:`[\x,\p]=i\hbar` Returns: tuple: the mean and the covariance matrix of the Gaussian state """# pylint: disable=unused-argument# Note: the internal order of mu and cov is different to the one used in Strawberry Fieldsreturncov,mu
[docs]defset_state(state,wire,cov,mu):r"""Inserts a single mode Gaussian into the state representation of the complete system. Args: state (tuple): contains covariance matrix and means vector of existing state wire (Wires): wire corresponding to the new Gaussian state cov (array): covariance matrix to insert mu (array): vector of means to insert Returns: tuple: contains the vector of means and covariance matrix. """cov0=state[0]mu0=state[1]N=len(mu0)//2# insert the new state into the means vectormu0[[wire[0],wire[0]+N]]=mu# insert the new state into the covariance matrixind=np.concatenate([wire.toarray(),wire.toarray()+N])rows=ind.reshape(-1,1)cols=ind.reshape(1,-1)cov0[rows,cols]=covreturncov0,mu0
[docs]defphoton_number(cov,mu,params,hbar=2.0):r"""Calculates the mean photon number for a given one-mode state. Args: cov (array): :math:`2\times 2` covariance matrix mu (array): length-2 vector of means params (None): no parameters are used for this expectation value hbar (float): (default 2) the value of :math:`\hbar` in the commutation relation :math:`[\x,\p]=i\hbar` Returns: tuple: contains the photon number expectation and variance """# pylint: disable=unused-argumentex=(np.trace(cov)+mu.T@mu)/(2*hbar)-1/2var=(np.trace(cov@cov)+2*mu.T@cov@mu)/(2*hbar**2)-1/4returnex,var
[docs]defhomodyne(phi=None):"""Function factory that returns the Homodyne expectation of a one mode state. Args: phi (float): the default phase space axis to perform the Homodyne measurement Returns: function: A function that accepts a single mode means vector, covariance matrix, and phase space angle phi, and returns the quadrature expectation value and variance. """ifphiisnotNone:def_homodyne(cov,mu,params,hbar=2.0):"""Arbitrary angle homodyne expectation."""# pylint: disable=unused-argumentrot=rotation(phi)muphi=rot.T@mucovphi=rot.T@cov@rotreturnmuphi[0],covphi[0,0]return_homodynedef_homodyne(cov,mu,params,hbar=2.0):"""Arbitrary angle homodyne expectation."""# pylint: disable=unused-argumentrot=rotation(params[0])muphi=rot.T@mucovphi=rot.T@cov@rotreturnmuphi[0],covphi[0,0]return_homodyne
[docs]defpoly_quad_expectations(cov,mu,wires,device_wires,params,hbar=2.0):r"""Calculates the expectation and variance for an arbitrary polynomial of quadrature operators. Args: cov (array): covariance matrix mu (array): vector of means wires (Wires): wires to calculate the expectation for device_wires (Wires): corresponding wires on the device params (array): a :math:`(2N+1)\times (2N+1)` array containing the linear and quadratic coefficients of the quadrature operators :math:`(\I, \x_0, \p_0, \x_1, \p_1,\dots)` hbar (float): (default 2) the value of :math:`\hbar` in the commutation relation :math:`[\x,\p]=i\hbar` Returns: tuple: the mean and variance of the quadrature-polynomial observable """Q=params[0]# HACK, we need access to the Poly instance in order to expand the matrix!# TODO: maybe we should make heisenberg_obs a class method or a static method to avoid this being a 'hack'?op=qml.PolyXP(Q,wires=wires)Q=op.heisenberg_obs(device_wires)ifQ.ndim==1:d=np.r_[Q[1::2],Q[2::2]]returnd.T@mu+Q[0],d.T@cov@d# convert to the (I, x1,x2,..., p1,p2...) orderingM=np.vstack((Q[0:1,:],Q[1::2,:],Q[2::2,:]))M=np.hstack((M[:,0:1],M[:,1::2],M[:,2::2]))d1=M[1:,0]d2=M[0,1:]A=M[1:,1:]d=d1+d2k=M[0,0]d2=2*A@mu+dk2=mu.T@A@mu+mu.T@d+kex=np.trace(A@cov)+k2var=2*np.trace(A@cov@A@cov)+d2.T@cov@d2modes=np.arange(2*len(device_wires)).reshape(2,-1).Tgroenewald_correction=np.sum([np.linalg.det(hbar*A[:,m][n])forminmodesforninmodes])var-=groenewald_correctionreturnex,var
[docs]deffock_expectation(cov,mu,params,hbar=2.0):r"""Calculates the expectation and variance of a Fock state probability. Args: cov (array): :math:`2N\times 2N` covariance matrix mu (array): length-:math:`2N` vector of means params (Sequence[int]): the Fock state to return the expectation value for hbar (float): (default 2) the value of :math:`\hbar` in the commutation relation :math:`[\x,\p]=i\hbar` Returns: tuple: the Fock state expectation and variance """# pylint: disable=unused-argumentex=fock_prob(cov,mu,params[0],hbar=hbar)# var[|n><n|] = E[|n><n|^2] - E[|n><n|]^2 = E[|n><n|] - E[|n><n|]^2var=ex-ex**2returnex,var
[docs]defidentity(*_,**__):r"""Returns 1. Returns: tuple: the Fock state expectation and variance """return1,0
[docs]classDefaultGaussian(Device):r"""Default Gaussian device for PennyLane. Args: wires (int, Iterable[Number, str]): Number of subsystems represented by the device, or iterable that contains unique labels for the subsystems as numbers (i.e., ``[-1, 0, 2]``) or strings (``['ancilla', 'q1', 'q2']``). Default 1 if not specified. shots (None, int): How many times the circuit should be evaluated (or sampled) to estimate the expectation values. If ``None``, the results are analytically computed and hence deterministic. hbar (float): (default 2) the value of :math:`\hbar` in the commutation relation :math:`[\x,\p]=i\hbar` """name="Default Gaussian PennyLane plugin"short_name="default.gaussian"pennylane_requires=__version__version=__version__author="Xanadu Inc."_operation_map={"Identity":Identity.identity_op,"Snapshot":None,"Beamsplitter":beamsplitter,"ControlledAddition":controlled_addition,"ControlledPhase":controlled_phase,"Displacement":displacement,"QuadraticPhase":quadratic_phase,"Rotation":rotation,"Squeezing":squeezing,"TwoModeSqueezing":two_mode_squeezing,"CoherentState":coherent_state,"DisplacedSqueezedState":displaced_squeezed_state,"SqueezedState":squeezed_state,"ThermalState":thermal_state,"GaussianState":gaussian_state,"InterferometerUnitary":interferometer_unitary,}_observable_map={"NumberOperator":photon_number,"QuadX":homodyne(0),"QuadP":homodyne(np.pi/2),"QuadOperator":homodyne(None),"PolyXP":poly_quad_expectations,"FockStateProjector":fock_expectation,"Identity":identity,}_circuits={}def__init__(self,wires,*,shots=None,hbar=2,analytic=None):super().__init__(wires,shots,analytic=analytic)self.eng=Noneself.hbar=hbarself._debugger=Noneself.reset()
[docs]defapply(self,operation,wires,par):# translate to wire labels used by devicedevice_wires=self.map_wires(wires)ifoperation=="Displacement":self._state=displacement(self._state,device_wires.labels[0],par[0]*cmath.exp(1j*par[1]))return# we are done hereifoperation=="GaussianState":iflen(device_wires)!=self.num_wires:raiseValueError("GaussianState covariance matrix or means vector is ""the incorrect size for the number of subsystems.")self._state=self._operation_map[operation](*par,hbar=self.hbar)return# we are done hereifoperation=="Snapshot":ifself._debuggerandself._debugger.active:gaussian={"cov_matrix":self._state[0].copy(),"means":self._state[1].copy()}self._debugger.snapshots[len(self._debugger.snapshots)]=gaussianreturn# we are done hereif"State"inoperation:# set the new device statecov,mu=self._operation_map[operation](*par,hbar=self.hbar)# state preparations only act on at most 1 subsystemself._state=set_state(self._state,device_wires[:1],cov,mu)return# we are done here# get the symplectic matrixS=self._operation_map[operation](*par)# expand the symplectic to act on the proper subsystemS=self.expand(S,device_wires)# apply symplectic matrix to the means vectormeans=S@self._state[1]# apply symplectic matrix to the covariance matrixcov=S@self._state[0]@S.Tself._state=[cov,means]
[docs]defexpand(self,S,wires):r"""Expands a Symplectic matrix S to act on the entire subsystem. Args: S (array): a :math:`2M\times 2M` Symplectic matrix wires (Wires): wires of the modes that S acts on Returns: array: the resulting :math:`2N\times 2N` Symplectic matrix """ifself.num_wires==1:# total number of wires is 1, simply return the matrixreturnSN=self.num_wiresw=wires.toarray()M=len(S)//2S2=np.identity(2*N)ifM!=len(wires):raiseValueError("Incorrect number of subsystems for provided operation.")S2[w.reshape(-1,1),w.reshape(1,-1)]=S[:M,:M].copy()# XXS2[(w+N).reshape(-1,1),(w+N).reshape(1,-1)]=S[M:,M:].copy()# PPS2[w.reshape(-1,1),(w+N).reshape(1,-1)]=S[:M,M:].copy()# XPS2[(w+N).reshape(-1,1),w.reshape(1,-1)]=S[M:,:M].copy()# PXreturnS2
[docs]defexpval(self,observable,wires,par):ifobservable=="PolyXP":cov,mu=self._stateev,var=self._observable_map[observable](cov,mu,wires,self.wires,par,hbar=self.hbar)else:cov,mu=self.reduced_state(wires)ev,var=self._observable_map[observable](cov,mu,par,hbar=self.hbar)ifself.shotsisnotNone:# estimate the ev# use central limit theorem, sample normal distribution once, only ok if n_eval is large# (see https://en.wikipedia.org/wiki/Berry%E2%80%93Esseen_theorem)ev=np.random.normal(ev,math.sqrt(var/self.shots))returnev
[docs]defsample(self,observable,wires,par):"""Return a sample of an observable. .. note:: The ``default.gaussian`` plugin only supports sampling from :class:`~.X`, :class:`~.P`, and :class:`~.QuadOperator` observables. Args: observable (str): name of the observable wires (Wires): wires the observable is to be measured on par (tuple): parameters for the observable Returns: array[float]: samples in an array of dimension ``(n, num_wires)`` """iflen(wires)!=1:raiseValueError("Only one mode can be measured in homodyne.")ifobservable=="QuadX":phi=0.0elifobservable=="QuadP":phi=np.pi/2elifobservable=="QuadOperator":phi=par[0]else:raiseNotImplementedError(f"default.gaussian does not support sampling {observable}")cov,mu=self.reduced_state(wires)rot=rotation(phi)muphi=rot.T@mucovphi=rot.T@cov@rotstdphi=math.sqrt(covphi[0,0])meanphi=muphi[0]returnnp.random.normal(meanphi,stdphi,self.shots)
[docs]defreset(self):"""Reset the device"""# init the state vector to |00..0>self._state=vacuum_state(self.num_wires,self.hbar)
[docs]defreduced_state(self,wires):r"""Returns the covariance matrix and the vector of means of the specified wires. Args: wires (Wires): requested wires Returns: tuple (cov, means): cov is a square array containing the covariance matrix, and means is an array containing the vector of means """iflen(wires)==self.num_wires:# reduced state is full statereturnself._state# translate to wire labels used by devicedevice_wires=self.map_wires(wires)# reduce rho down to specified subsystemsind=np.concatenate([device_wires.toarray(),device_wires.toarray()+self.num_wires])rows=ind.reshape(-1,1)cols=ind.reshape(1,-1)returnself._state[0][rows,cols],self._state[1][ind]
[docs]defexecute(self,operations,observables):iflen(observables)>1:raiseqml.QuantumFunctionError("Default gaussian only support single measurements.")returnsuper().execute(operations,observables)