Exact Boson Sampling using Gaussian continuous variable measurements

TL;DR

This paper introduces a method for exact Boson Sampling using Gaussian continuous-variable measurements, demonstrating classical hardness for exact sampling, but not for approximate sampling, thus advancing quantum computational complexity understanding.

Contribution

It presents a novel device setup combining Fock states, linear interactions, and Gaussian measurements that achieves classically hard exact sampling, extending Boson Sampling concepts.

Findings

Exact sampling is classically hard for the proposed device.

Gaussian measurements do not imply hardness for approximate sampling.

The paper discusses conditions needed for approximate sampling hardness.

Abstract

BosonSampling is a quantum mechanical task involving Fock basis state preparation and detection and evolution using only linear interactions. A classical algorithm for producing samples from this quantum task cannot be efficient unless the polynomial hierarchy of complexity classes collapses, a situation believe to be highly implausible. We present method for constructing a device which uses Fock state preparations, linear interactions and Gaussian continuous-variable measurements for which one can show exact sampling would be hard for a classical algorithm in the same way as Boson Sampling. The detection events used from this arrangement does not allow a similar conclusion for the classical hardness of approximate sampling to be drawn. We discuss the details of this result outlining some specific properties that approximate sampling hardness requires.

Figures

Equations

$$F=\sum _{x\in Q}f(x)$$

$$F{g}^{-1}\le \tilde{F}\le Fg$$

$$F-g\le \tilde{F}\le F+g.$$

$${\Pi }^n(x_{1,\theta },p_{2,\theta })=\frac{1}{2\pi }D(\alpha )|n⟩⟨n|D^{†}(\alpha ).$$

$${\Pi }^n(R)=n!e^{-R}\sum _{k=0}^{\infty }\frac{R^{k-n}}{k!}{\left(L_n^{k-n}(R)\right)}^2|k⟩⟨k|$$

$${\Pi }_t^1$$

$$=\sum _{k=0}^{\infty }G(t,k)|k⟩⟨k|,$$

$$G(t,k)=\frac{1}{k!}\left(k^2\gamma (k,t)-2k\gamma (1+k,t)+\gamma (2+k,t)\right)$$

$${\Pi }_{\bar{t}}^1$$

$${\Pi }_t^1=\left(\frac{t^2}{2}-\frac{t^3}{3}+O(t^4)\right)|0⟩⟨0|$$

$$+\left(t-\frac{3t^2}{2}+\frac{7t^3}{6}+O(t^4)\right)|1⟩⟨1|+\left(t^2-\frac{4t^3}{3}+O(t^4)\right)|2⟩⟨2|$$

$$+\left(\frac{1}{2}{t}^3+O(t^4)\right)|3⟩⟨3|+O(t^4)|4⟩⟨4|+\cdots$$

$${\Pi }_{ϵ}^1=ϵ|1⟩⟨1|.$$

$$\eta (t)$$

$$p_D(t)$$

$$|1_N⟩=|\underset{N}{\underbrace{1,1,\cdots ,1}},\underset{M-N}{\underbrace{0,0,\cdots ,0}}⟩.$$

$${\Pi }^1(\mathbf{\alpha })=\frac{1}{(2\pi {)}^M}D(\mathbf{\alpha })|1_M⟩⟨1_M|D^{†}(\mathbf{\alpha }),$$

$$P(\mathbf{\alpha })$$

$$=\frac{1}{(2\pi {)}^M}{\left|\sum _{\mathbf{n}}⟨1_M|D^{†}(\mathbf{\alpha })|\mathbf{n}⟩⟨\mathbf{n}|{\mathcal{U}}_{\text{LON}}|1_N⟩\right|}^2.$$

$$⟨\mathbf{n}|{\mathcal{U}}_{\text{LON}}|1_N⟩=\text{Per}\left(U_{1_N\times \mathbf{n}}\right),$$

$$b_j^{†}={\mathcal{U}}_{\text{LON}}{a}_j^{†}{\mathcal{U}}_{\text{LON}}^{†}=\sum _{i=1}^M{U}_{ij}{a}_i^{†},$$

$$⟨1_M|D^{†}(\mathbf{\alpha })|\mathbf{n}⟩=\prod _{j=1}^M{e}^{-|{\alpha }_j{|}^2/2}\frac{({\alpha }_j^{\ast }{)}^{n_j-1}}{\sqrt{n_j!}}(n_j-|{\alpha }_j{|}^2).$$

$$P(\mathbf{\alpha })=\frac{e^{-\mathbf{\alpha }{\mathbf{\alpha }}^{†}}}{(2\pi {)}^M}{\left|\sum _{\mathbf{n}}\text{Per}\left(U_{1_N\times \mathbf{n}}\right)\prod _{j=1}^M\frac{({\alpha }_j^{\ast }{)}^{n_j-1}}{\sqrt{n_j!}}(n_j-|{\alpha }_j{|}^2)\right|}^2.$$

$$P({\mathbf{\alpha }}_{{\mathbf{n}}_N})=\frac{1}{(2\pi {)}^M}{\left|\text{Per}\left(U_{1_N\times {\mathbf{n}}_N}\right)\right|}^2\prod _{n_{N,j}=0}{e}^{-|{\alpha }_j{|}^2}|{\alpha }_j{|}^2,$$

$$P({\mathbf{\alpha }}_{1_N})$$

$$=\frac{1}{(2\pi {)}^M}{\left|\text{Per}\left(U_{1_N\times 1_N}\right)\right|}^2\prod _{j=N+1}^M{e}^{-|{\alpha }_j{|}^2}|{\alpha }_j{|}^2.$$

$$P_{\text{P}}(\mathbf{R})=e^{-|\mathbf{R}|}$$

$$P_{\text{P}}({\mathbf{R}}_{{\mathbf{n}}_N})=$$

$$\times \prod _{n_{Nj}=0}\frac{1}{n_j!}{e}^{-R_j}{R}_j^{n_j-1}(n_j-R_j{)}^2.$$

$$P_{\text{D}}({\mathbf{m}}_N)=$$

$$\times \prod _{m_{Nj}=1}G(t,n_j)\prod _{m_{Nj}=0}(1-G(t,n_j)).$$

$$P_{\text{D}}({\mathbf{m}}_N)=$$

$$+\sum _{\mathbf{m}={\mathbf{m}}_N^1}|\text{Per}\left(U_{1_N\times \mathbf{m}}\right){|}^{2}O(t^{N+2}),$$

$$P_{\text{D}}(1_N)=|\text{Per}\left(U_{1_N\times 1_N}\right){|}^2{t}^N.$$

Full text

I Introduction

$\BosonSampling$ is the task of producing statistical samples from Fock basis measurements of a bosonic $M$-mode linear scattering network with an input consisting of $N$ modes prepared with a single boson and the remaining $M-N$ modes prepared in the vacuum state. This task, whilst not universal for quantum computing, has been shown to not be efficiently computable by any classical algorithm (or the polynomial hierarchy of complexity classes collapses which is believe to be extremely unlikely) AA . However, a quantum implementation is efficient as one merely needs to build the scattering device as described within a sufficiently small error budget. This is many orders of magnitude easier than the construction of a fully universal quantum computer, but is still a challenge for current technology.

Attempts have been made to identify scenarios where the proof of classical hardness of $\BosonSampling$ can be used or adapted to other sampling problems. One particular scenario that is of experimental interest is in the use of continuous variable (CV) Gaussian states or measurements. For another restricted computational model based on sampling from qubit circuits involving commuting coherent rotations it has been shown that CV variants are hard to simulate classically Douce17 . It is known that for linear networks with Gaussian state inputs and Gaussian measurements it is efficient for classical algorithms to not only produce samples but compute the entire output distribution Bartlett03 . Never-the-less it has been shown that a hybrid approach which involves linear networks with input two-mode squeezed vacuum states and Fock basis detection has a similar classical hardness proof to the original $\BosonSampling$ problem Lund2014 . There is also evidence for the classical hardness of a more general construction involving squeezed-vacuum inputs, linear optics and Fock basis detection Rahimi2015 ; Ham2016 . The question this paper addresses is the reverse situation, Fock state inputs to linear networks and homodyne detection.

An important aspect of the hardness proof for $\BosonSampling$ is that the output probability distribution contains probabilities which are proportional to matrix permanents from sub-matricies of the matrix describing the linear scattering network. A matrix permanent is a quantity which is computed like a matrix determinant without the alternating addition and subtraction. In fact, when sampling with Fock state inputs and detection, all detection probabilities are proportional to sub-matricies derived from the linear scattering network Scheel04 . This special situation, given two plausible conjectures hold, allows for the hardness proof of “approximate” sampling to be shown AA . This is because an allowed error budget’s effect can be spread over all detection events provided the linear network appears randomly distributed.

Here we show that using single-photon input states and a particular Gaussian measurement one can extract sub-matrix permanents to within an exponentially small error. Therefore one can show that exact sampling from this distribution is hard. It is not necessarily the case that approximate sampling is still hard and we discuss this in relation to our construction.

In section II we present some of the background behind the hardness arguments for $\BosonSampling$. Then in section III we present the CV-n detector model using a Fock state $\ket{n}$ input with CV measurements. In section IV we will then describe exact sampling using the CV-1 model and the technical details involved in showing the hardness of computing samples from the CV output distribution. Finally we will discuss the issues preventing the hardness result from being used in this model to make definitive conclusions about the hardness or not of approximate sampling.

II Classical hardness of $\BosonSampling$

A problem in the class $\BosonSampling$ is one where the statistical samples can be generated by an $M$ mode linear interaction between $N$ singly occupied bosonic modes (and $M-N$ bosonic vacua) which is subsequently detected in the Fock basis. $\BosonSampling$ is either inefficient using classical computational resources (not in $\P$), or the “polynomial hierarchy” of complexity classes collapses to the third level, a situation believe to be implausible. It is not our goal to present in full the background and subsequent arguments towards the truth of this statement as this has been done elsewhere AA . We will, however, outline some of the key aspects used in this paper that are needed to understanding what is required for the proof presented in AA to hold.

II.1 Polynomial hierarchy

The polynomial hierarchy of complexity classes is a nested structure defined by the use of oracles. An oracle is essentially an assumption on the resources available to an algorithm which can greatly assist in proving statements in computational complexity. The hierarchy has a complex definition and we will concentrate on a simplified version.

The class $\NP$ is the set of decision (’yes’/’no’) problems whose satisfying input can be verified efficiently. This defines the first level of the polynomial hierarchy. The second level is then the class $\NP$ with access to an oracle from the first level. Subsequent levels are defined by continuing this recursion, e.g. the third level is the class $\NP$ with access to an oracle from the second level.

This structure has a strong connection to similarly defined hierarchies within number theory and set theory. In those cases each level of the hierarchy is strictly larger than lower levels. If two levels were to coincide, then the addition of levels stops growing the hierarchy is said to collapse. In terms of the computational complexity structure, a collapse of the hierarchy to the first level means that $\P=\NP$ or that $\NP$ problems can be efficiently solved deterministically a situation believed to be highly implausible. A collapse to the second level would mean $\P^{\NP}=\NP^{\NP}$ which is the same statement relative to an $\NP$ oracle. Being relative to the oracle means that the statement is slightly more plausible, but it is still believed to be not possible. The prevailing belief is that the polynomial hierarchy does not collapse to any level and this is the assumption on which the hardness of $\BosonSampling$ can be proven.

II.2 Stockmeyer’s approximate counting algorithm

Critical to the hardness proof of $\BosonSampling$ is the use of Stockmeyer’s approximate counting algorithm Stockmeyer1983 . This algorithm computes estimates of a quantity defined as

[TABLE]

where $f:Q\rightarrow\{0,1\}$ is a boolean function from length $l$ bitstrings $Q=\{0,1\}^{l}$ onto a single bit. In other words $F$ is the number of inputs that result in an output of $1$, a set we will call $Q_{1}$. The computed estimate of this quantity is multiplicative, which means the estimate of $F$ that the algorithm produces is $\tilde{F}$ satisfying

[TABLE]

where $g>1$ and for Stockmeyer’s algorithm is lower bounded by $1+1/poly(l)$.

Stockmeyer’s algorithm computes the estimate $\hat{F}$ by finding the smallest output with no collisions for a randomly chosen hash functions on $Q_{1}$. That is, choose randomly a function $h:Q_{1}\rightarrow\{0,1\}^{p}$ (for $p\leq l$) and if there exists no elements $a,b\in Q_{1}$ such that $h(a)=h(b)$ (a hash collision) then we have made an upper bounded estimate on the size of $Q_{1}$ of $2^{p}$.

Here we can see the elements which make up this algorithm: the estimate of size is multiplicative, it requires finding hash collisions, which is an $\NP$ problem, and involves random choices. This results in Stockmeyer’s algorithm being contained within a class $\FBPP^{\NP}$ (if the function $f$ is efficiently computable). The superscript notation describes access to an $\NP$ oracle, which in the case of Sockmeyer’s algorithm is used to find hash collisions. $\BPP$ means the algorithm proceeds by making random (probabilistic) choices with a probability of success at least $2/3$. The prefix ’F’ is used to describe the output from this algorithm, which is a number (function) rather than a decision or ’yes’/’no’ output. This is therefore not an algorithm one would expect to be efficiently computable when realistically implemented. But it is used here to establish if a problem lives within the “polynomial” hierarchy of complexity classes. In this case, this algorithm lives within the third level as $\BPP\subset\NP^{\NP}$ Sip ; Laut .

This algorithm is used to make multiplicative estimations of the underlying probabilities of a distribution from samples of the distribution. In the model of classical computation, unlike quantum mechanical models, there is no inherent randomness. Randomness is introduced by external means and can be regarded as an input to the algorithm. If the function $f(x)$ above represent a sampling algorithm for a two-outcome distribution, then it can be thought of as converting input random bits $x$ distributed uniformly into samples of the desired probability distribution. One can then estimate the probability of the outcome $1$ by counting how many inputs produce the outcome $1$ relative to the total number of possible inputs. This would mean dividing the estimate $\hat{F}$ by $2^{l}$ which will also produce a multiplicative estimate even though it has been divided by an exponentially increasing factor.

A polynomial hierarchy collapse is triggered if $f$ is efficient to compute and the probability it produces samples from is a quantity for which multiplicative estimation is known to be outside the third level of the polynomial hierarchy.

II.3 Multiplicative and additive errors

Producing approximations within multiplicative errors is particularly powerful. It is more natural to consider the case of additive errors. An additive estimate of $F$ would be one satisfying

[TABLE]

This kind of estimation will arise from more natural models of errors within a computation. For example, the function $f$ above most likely admits errors which are of an additive nature and hence the estimations performed using this noisy model will also be additive. Only in the case of admitting exactly zero noise to $f$, or in the sampling case a situation called “exact sampling”, can one generate a multiplicative estimation.

The crucial outcome of AA is that the polynomial hierarchy collapse can be shown to occur in $\BosonSampling$ even with a given level of total variation distance between the ideal and actual distributions. The total variation distance is an additive quantity and will generate additive errors in estimates. But as the probabilities in $\BosonSampling$ tend to decrease to zero exponentially, the introduction of additive errors will overwhelm the magnitude of the quantity being estimated. The trick of AA is to use the structure of the probabilities in $\BosonSampling$ to encode the estimated quantity within an exponentially large set of possible outcomes in a way that without knowledge of where the problem is encoded, looks like a Haar-random unitary matrix (i.e. the choice is hidden from the implementation). This means that the additive estimate generated will have bounds determined by the average case which turns out to be sufficient to additively estimate matrix permanents. Using this estimate it is then possible to trigger a polynomial hierarchy collapse to the third level, given that two plausible conjectures about the nature of estimating permanents of Gaussian random matrices hold true.

III The detector

We will now describe our continuous-variable detection model and how we use it to construct a probability distribution which can be used in the arguments of AA . The model is based on a measurement device for measuring in the displaced number state basis, which we refer to it as CV-n measurement, and two variations of this measurement, phase-randomised CV-n (PRCV-n) and discretised-phase-randomised CV-n (DPRCV-n). Only a brief outline of the model is presented here. The details of this calculation are given in Appendix A.

As depicted in Fig. 1, the measurement device works as follows. The input signal is overlapped on a 50:50 beamsplitter with a number state $\ket{n}$, and the outputs of the beamsplitter are measured by two conjugate homodyne detections whose local oscillators have $\pi/2$ phase difference. As shown in Appendix A, the POVM elements of this measurement are

[TABLE]

The two real numbers $x_{1,\theta},p_{2,\theta}$, form the results of a simultaneous measurement of two orthogonal quadratures. Notice that for $n=0$, CV-0, we have heterodyne measurement.

Phase-randomised CV-n (PRCV-n). If the phase $\theta$ of the local oscillators is randomised, while the relative phase is fixed, we have PRCV-n measurement. As discussed in Appendix A, the outcome of this measurement is a single non-negative parameter $R=r^{2}=x_{1,\theta}^{2}+p_{2,\theta}^{2}$ with the corresponding POVM element

[TABLE]

with $L^{m}_{k}(.)$ being the generalized Laguerre polynomials. Though this measurement has less information, we will concentrate on it as we are only interested in events where $R\approx 0$ and having the POVM diagonal in the Fock basis greatly simplifies many calculations.

Discretized-phase-randomised CV-1 (DPRCV-1). The measurement model outputs continuous variables (i.e. $R$) and hence these POVMs are describing probability densities. To be able to use the methods of AA to make some definitive statement about computational hardness, one must work in probabilities and not probability densities. Hence, we will utilize the discretized version of the phase-randomized CV-1 (DPRCV-1) measurement through its ability to distinguish between the single-photon state and other number states. We divide the range of $0\leq R\leq\infty$ for the POVM elements in Eq. (5) for $n=1$, to two parts $t=\{R|0\leq R\leq t\}$ which represents the smallest discrete region and $\bar{t}=\{R|t\leq R\leq\infty\}$ which represents all the others. Thus, we have two POVM elements for each interval:

[TABLE]

where

[TABLE]

with $\gamma(k,t)=\int_{0}^{t}dRe^{-R}R^{k-1}$ being the lower incomplete Gamma function. As this is a two-outcome measurement and each of the elements is bounded above by $1$ the POVM for all other outcomes can then be written

[TABLE]

The POVM element in Eq. (6) can be expanded in a power series as

[TABLE]

Therefore, if $t$ is small $t=\epsilon$ such that $\epsilon^{2}\approx 0$, the POVM element associated with detecting $R\in[0,\epsilon]$ is

[TABLE]

It is this form of the POVM which allows for a measurement that distinguishes a single photon Fock state from the remainder of the Hilbert space.

To demonstrate what this measurement is detecting consider the case of an input states being either $\ket{0}$ or $\ket{1}$. Then when detecting the single photon state (and for any value of $t$),

[TABLE]

can be thought of as the efficiency of the detector. Also

[TABLE]

can be thought as the dark count probability. Fig. 2 compares these two quantities.

IV Boson sampling

We now consider the problem of sampling from the output probability distribution of a linear-optical network (LON) using CV-1, PRCV-1, and DPRCV-1 measurements (where the ancillary Fock state is $n=1$) introduced in the previous section. We first derive the output probability distribution for each measurement scheme, and then discuss the complexity of sampling from the probability distribution.

IV.1 Probability distributions

In this setup, the input state to the LON is

[TABLE]

Each output mode is measured by a CV-1 measurement; hence, the overall POVM elements are

[TABLE]

where $\bm{\alpha}=(\alpha_{1},\alpha_{2},\cdots,\alpha_{M})$. As shown in Fig. 3, this setup is equivalent to injecting single photons into a larger network and performing homodyne measurements at the output, i.e., boson sampling with homodyne measurements.

Therefore, the output probability density is given by

[TABLE]

Here $\mathcal{U}_{\text{LON}}$ is the unitary operation associated with the LON, and $\ket{\bm{n}}=\ket{n_{1},n_{2},\cdots,n_{M}}$ is the multimode Fock state. As a LON preserves the number of photons, the sum is restricted to $\bm{n}$’s which satisfy $\sum_{i=1}^{M}n_{i}=N$.

We have Scheel04

[TABLE]

which is the permanent of an $N\times N$ submatrix of the unitary matrix of LON, $U_{1_{N}\times\bm{n}}$ that is corresponding to the first $N$ rows and columns with multiplicity of $n_{i}$. The unitary matrix $U$ is defined through the relation

[TABLE]

where $b^{\dagger}_{j}$ and $a^{\dagger}_{i}$ are modal creation operators for output $j$ and input $i$, respectively. Also, using the expression for $\bra{k}D(\alpha)\ket{m}$ in Cah-Glau , one can simply verify that

[TABLE]

By using these relations, the output probability density (15) becomes

[TABLE]

From this expression we can see that if $\alpha_{j}=0$, the probability density $P(\bm{\alpha})$ is zero unless $n_{j}=1$.

We define $\bm{n}_{N}=(n_{N,1},n_{N,2},\cdots,n_{N,M})$ as a $M$-tuple whose elements are either zero or one and the number of ones is $N$, and define $\bm{\alpha}_{\bm{n}_{N}}$ a vector of length $M$ whose element are given by $[\bm{\alpha}_{\bm{n}_{N}}]_{j}=(1-n_{N,j})\alpha_{j}$, i.e., it has $N$ zero elements corresponding to $n_{N,j}=1$. Using these notations, we can see that

[TABLE]

where the product is over $M-N$ nonzero elements of $\bm{\alpha}_{\bm{n}_{N}}$. This implies that if $N$ of $M$ $\alpha$’s are zero, the probability density $P(\bm{\alpha})$ is proportional to the absolute value square of permanent of a submatix of $U$. For example, for $\bm{n}_{N}=1_{N}$

[TABLE]

If we use PRCV-1 measurements at the output of LON, using the POVM (68), the output probability density is given by

[TABLE]

where $\bm{R}=(R_{1},R_{2},\cdots,R_{M})$ and $|\bm{R}|=R_{1}+R_{2}+\cdots+R_{M}$. Defining $\bm{R}_{\bm{n}_{N}}$ similar to $\bm{\alpha}_{\bm{n}_{N}}$, we have

[TABLE]

For DPRCV-1 measurements at the output, the probability of detecting outcome $\bm{m}_{N}$, that is, $N$ clicks within interval $t$ and $M-N$ clicks within $\bar{t}$, using Eq. (7), is given by

[TABLE]

If $t$ is very small, this probability can be expanded to leading order in a powers of $t$ as

[TABLE]

where $\bm{m}_{N}^{1}$ are $M$-tuples whose elements are equal to $\bm{m}_{N}$ except one 0 and one 1 are interchanged.

Notice that the probability distribution (23) corresponding to DPRCV-1 measurements is coarse-grained version of probability density (21), and that itself is phase randomised version of (18). Therefore, if exact sampling from the probability distribution (23) is classically hard, exact sampling from the other probability distributions must be classically hard as well. In the following two subsections we discuss exact and approximate sampling from the probability distribution (23).

IV.2 Exact boson sampling

To leading order in $t$, Eq. (24) give a distribution for which hardness of exact sampling can be determined as we shall show in this section. Merely setting the $O(t^{N+1})$ terms to zero results in a probability

[TABLE]

A probability of this form does not pose any problems for the argument of hardness in exact sampling. However, it is important to consider higher order terms as exact sampling allows for sampling using any matrix even those with vanishingly small permanents. For example, one could consider the case where the permanent $|\text{Per}\left(U_{1_{N}\times 1_{N}}\right)|^{2}=O(t^{2})$ but the unwanted higher order terms have $|\text{Per}\left(U_{1_{N}\times\bm{n}}\right)|^{2}=O(1)$. Furthermore, as $|\text{Per}\left(U_{1_{N}\times\bm{n}}\right)|^{2}$ defines a probability distribution over all $\bm{n}$ (i.e. the $\BosonSampling$ distribution) we have

[TABLE]

This means that, using Eq. (24), the probability can be written as

[TABLE]

where $E=O(t^{2})$ is an “error” over the desired probability. Unfortunately $E$ contributes to the probability in an additive sense. However, as we have used the upper bound on the matrix permanents that contribute to $E$ we know that this bound is independent of $N$. So if we choose $E$ to be a small constant, we need to quadratically vary $t$ to achieve that constant. Furthermore, with constant $E$, the sizes of sub-matrix permanents that can be estimated must be lower bounded with a bound that depends on $E$.

To make this more explicit consider the case where $0<L\leq|\text{Per}\left(U_{1_{N}\times 1_{N}}\right)|^{2}\leq 1$. If we used the approximate counting algorithm using $P_{\text{D}}(1_{N})$ from Eq. (27) then we would have an estimate $\tilde{p}$ which satisfies

[TABLE]

where $g>1$ is the multiplicative error factor from the estimation. We can also write

[TABLE]

and using $|\text{Per}\left(U_{1_{N}\times 1_{N}}\right)|^{2}>L$

[TABLE]

Now if $|E|/L<1/2$ (arbitrarily, could be any constant less than 1), we have

[TABLE]

where we now have a multiplicative error estimate with error $g^{\prime}=(1+2|E|/L)g$.

The #P-hardness of approximating $\text{Per}(X)^{2}$ (Theorem 28 from AA ) where $X$ is an $N\times N$ matrix, requires a choice of $g$ that is polynomially dependent on $N$. This will be achieved here if the extra $(1+2|E|/L)$ term maintains a polynomial scaling.

At this point, we consider the origin of the added error term $E$ to argue how $L$ must scale in terms of $N$. All of the schemes described in Section III depend on a continuous output parameter based on the CV-1 style measurement. However, any experimental realisation will output values to within some precision. One particular way to do this, with a connection to discrete computational outputs is to consider a discretisation of the results into $b$-bit integers. The outcomes $R$ from PRCV-1 are non-negative, so in this case consider $b$-bit non-negative integers. $R$ is unbounded for large values, but in our scheme we are only interested in small values so we can choose an arbitrary fixed boundary above which all values are assigned the same bitstring output. With these requirements, results can be partitioned into $2^{b}$ regions such that the largest result corresponds to all outcomes above some fixed value, say $1$. In the simplest case, if $b=1$ then we would have the one bit binary string with ‘[math]’ representing all results of the continuous value $R$ between [math] and $1$ and the binary string ‘$1$’ representing values greater than $1$. Increasing $b$ and continuing this division of the results, we find that the range of values covered by equal size partitions of the values between [math] and $1$ we have $t=2/(2^{b}-1)=O(2^{-b})$. It is the scaling in $b$, and not $t$, that is required to be considered when analysing the complexity of the problems modelled on this device.

We will proceed with the analysis below assuming this simple partitioning of the values of $R$. It should be noted that other discretisation strategies can be considered. For example, a change of the discretisation thresholds such that smaller regions around zero can be considered by using a common technique called “companding”. However, one must be careful here not to sacrifice too much probability by having the zero region scale super-exponentially. Any scaling which results in the region at zero scaling as $O(2^{-poly(b)})$ is sufficient for the exact sampling argument below. This is due to the nature of the Stockmeyer approximate counting algorithm utilized in AA. As these strategies do not change the complexity hardness result we present, we will not consider them further in this paper.

Now we will return to considering the error term $|E|$. This quantity, as stated above, scales as $O(t^{2})$. But in the number of bits used for discretisation it will scale as as $O(2^{-2b})$. To achieve a $poly(N)$ scaling of $|E|/L$ a scaling in $L$ of $O(2^{-2b})$ would counteract the error leaving a $O(1)$ overhead.

So our procedure will achieve a multiplicative estimate with a very quickly decaying lower bound $L$ on the matrix permanents that can estimate. But what actually determines the choice of the lower bound $L$ that is permitted? The proof of the #P-hardness for approximating $\text{Per}(X)^{2}$ proceeds in AA by using this estimation polynomially many times to compute $\text{Per}(X)$ given that $X$ was a zero-one matrix. As part of the proof one needs to use the estimation procedure in a binary search towards a matrix whose permanent is zero and the procedure stops when sufficient precision has been achieved to determine the permanent of the zero-one matrix, as it must be an integer value. Introducing a lower bound on the permanents for which the estimation is valid would have to be compatible with this final precision.

Never-the-less, in the hardness proof for approximate $\BosonSampling$, it is assumed that the input matrices $X$ have matrix elements that have a Gaussian distribution. If one is to accept the conjectures of AA , then there is a low probability of randomly having $\text{Per}(X)$ near zero. Specifically, using the words from AA :

…if $X\sim\mathcal{G}^{n\times n}$ is Gaussian then …a $1-1/poly(n)$ fraction of [$\text{Per}(X)$’s] probability mass is greater than $\sqrt{n!}/poly(n)$ in absolute value, $\sqrt{n!}$ being the standard deviation.

Using this we could choose the lower bound of $Per(X)$ to scale as $\sqrt{N!}/poly(N)$. However, we have written the bound for our construction in terms of the submatrix of the unitary matrix $U$. Following the embedding procedure from AA this reduces the size of $X$ by a factor involving the matrix norm $||X||$. Therefore the probability reduces by a factor of $||X||^{2N}\leq 2^{poly(N)}$ as $N$ (the number of input photons) is a proxy for the problem input size (i.e. the matrix $X$) which must be efficiently represented in $N$. This results in a scaling of $L=O(||X||^{-2N}N!/poly(N))$. In terms of the discretisation bits the scaling is $b=O(N\log N+poly(N))$ where the $poly(N)$ in this scaling is from the polynomial which bounds $||X||^{-2N}$. So a polynomial scaling in the size of the discretisation allows for approximation of permanents which are highly probable when using an approximate $\BosonSampling$ algorithm. This argument shows that the exact sampling problem presented here, which permits multiplicative permanent estimation with permanents with an exponentially decreasing lower bound, must also be hard to compute with classical resources.

IV.3 Approximate sampling

The arguments just made do not permit one to continue the hardness argument through to the case of approximate sampling as was done in AA for the CV distributions we have constructed above. The reason is quite simple, only an exponentially small subset of events (i.e. detections near the origin) are used to make the exact sampling argument. Conversely, for the same reason, one also cannot definitively conclude that approximate sampling from this distribution is an efficient classical task.

The approximate sampling criterion for this CV distribution would be

[TABLE]

where $P$ is the probability density for our distribution of CV events from the device described above and $Q$ is the computable approximation. The event used for sampling above is the set of $\bm{\alpha}$, where $N$ $\alpha$’s are around a ball of radius $t$ around the origin. Hence all of this error could potentially be concentrated on our event. This would dominate any exponential pre-factors of the matrix permanent in the probability exponentially reducing the signal that is being fed into the approximate counting algorithm.

Approximate sampling is shown in AA by potentially utilising all possible (or more precisely the collision-free subspace of) events to perform the approximate counting. The algorithm cannot know which event is being probed, so a concentrate of error as described above would result in the approximate sampling algorithm being an exact sampling algorithm for almost all results.

So what is required in the CV case is more events. To show classical hardness for this distribution, one needs to argue that events away from the origin are also #P-hard when used exactly. This is difficult for this construction as these events do not reduce to Fock basis measurements under any approximation. Also, there must be combinatorially enough events so that the approximate sampling error can be considered low enough on average over all the events so that hardness is maintained. Even if this specific criterion cannot be met, this does not mean that the distribution is efficient for a classical computation. To show classical efficiency, either a constructive proof demonstrating efficiency (e.g. the methodology used in SRK16-suff ) or some other proof technique ruling out classical hardness would be required.

V Conclusion

We have constructed a continuous variable measurement model we call CV-n which measures in a displaced number state basis. The measurement model can be achieved by mixing a $n$ photon state with the input on a 50:50 beam splitter then making homodyne measurements on the output. We have shown how a $n$ Fock basis state measurement can be approximately achieved using this measurement utilising those cases when the homodyne measurement outcomes are simultaneously small.

We then showed that this measurement model is compatible with the exact $\BosonSampling$ problem, a computing task that has been shown to be inefficient for any classical device to simulate. We have discussed how this model as presented here is not compatible with approximate $\BosonSampling$ as the detection events utilised, when compared with the whole event space, is too small.

Acknowledgements

This work was supported by the Australian Research Council Centre of Excellence for Quantum Computation and Communications Technology (Project No. CE110001027).

Appendix A CV-n measurement: Continuous variable measurement in the displaced number state

In this appendix we obtain the POVM elements of the CV-n measurement; see Fig. 1. The POVM elements of this measurement are of this form

[TABLE]

where $x_{1,\theta}$ and $p_{2,\theta}$ are the outcomes of the fist and the second homodyne measurements, respectively, and $c$ is the normalization constant such that $\int dx_{1,\theta}\;dp_{2,\theta}\Pi(x_{1,\theta},p_{2,\theta})=I$.

The POVM elements of the first homodyne measurement are

[TABLE]

where

[TABLE]

Notice that $\ket{x_{1,\theta}=0}$ is an infinitely squeezed vacuum state, and can be written as

[TABLE]

with $\hat{S}(\xi)=\exp\left[(\xi^{*}\hat{a}^{2}-\xi\hat{a}^{\dagger 2})/2\right]$ being the squeezing operator and $\ket{0}$ being the vacuum state.

For the second homodyne measurement we have

[TABLE]

where

[TABLE]

Similarly, we have

[TABLE]

Now using Eqs. (34) and (37), $\ket{\Psi(x_{1,\theta},p_{2,\theta})}$ in Eq. (33) is given by

[TABLE]

where $\mathcal{N}$ is the normalization constant. Using the following relations for the 50:50 beamsplitter

[TABLE]

the expression $A$ in Eq. (40) becomes

[TABLE]

where

[TABLE]

and $D(\beta)$ is the displacement operator. Thus, Eq. (40) becomes

[TABLE]

Using

[TABLE]

we can write Eq. (46) as

[TABLE]

It can be shown that

[TABLE]

where

[TABLE]

and the second equality can be seen using this expression

[TABLE]

with $L^{m}_{k}(.)$ being the generalized Laguerre polynomials.

By substituting from Eq. (51) into Eq. (50), we get

[TABLE]

As this state must be normalized in the limit of $r\rightarrow\infty$, it is straightforward to see that

[TABLE]

As $\lim_{r\rightarrow\infty}\tanh r=1$, we now have

[TABLE]

where

[TABLE]

Therefore, POVM elements of the CV-n measurement are

[TABLE]

Notice that, it can be simply checked that

[TABLE]

A.1 Phase randomized CV-n measurement

If the phase $\theta$ of the local oscillator is randomized, with the relative phase $\pi/2$ being fixed, we have

[TABLE]

where

[TABLE]

The matrix elements of this operator in the Fock basis are

[TABLE]

Therefore, the POVM elements of the phase randomized measurement are given by

[TABLE]

It can be simply verified that

[TABLE]

For $n=1$, using

[TABLE]

we have

[TABLE]

so the POVM elements (68) becomes

[TABLE]

Notice that

[TABLE]

as expected.