The Classical Complexity of Boson Sampling

TL;DR

This paper presents a significantly faster classical algorithm for exact Boson Sampling, challenging the assumption that quantum supremacy in this domain is easily achievable and extending the feasible classical computation size.

Contribution

The authors introduce a new classical algorithm for exact Boson Sampling with exponential speedup over previous methods, enabling larger system sizes to be computed classically.

Findings

Classical algorithm runs in O(n 2^n + poly(m,n)) time

Algorithm can solve larger Boson Sampling instances than previous methods

Reduces the likelihood of near-term quantum supremacy in Boson Sampling

Abstract

We study the classical complexity of the exact Boson Sampling problem where the objective is to produce provably correct random samples from a particular quantum mechanical distribution. The computational framework was proposed by Aaronson and Arkhipov in 2011 as an attainable demonstration of `quantum supremacy', that is a practical quantum computing experiment able to produce output at a speed beyond the reach of classical (that is non-quantum) computer hardware. Since its introduction Boson Sampling has been the subject of intense international research in the world of quantum computing. On the face of it, the problem is challenging for classical computation. Aaronson and Arkhipov show that exact Boson Sampling is not efficiently solvable by a classical computer unless $P^{\#P} = BPP^{NP}$ and the polynomial hierarchy collapses to the third level. The fastest known exact classical algorithm for the standard Boson Sampling problem takes $O({m + n -1 \choose n} n 2^n )$ time to produce samples for a system with input size $n$ and $m$ output modes, making it infeasible for anything but the smallest values of $n$ and $m$. We give an algorithm that is much faster, running in $O(n 2^n + \operatorname{poly}(m,n))$ time and $O(m)$ additional space. The algorithm is simple to implement and has low constant factor overheads. As a consequence our classical algorithm is able to solve the exact Boson Sampling problem for system sizes far beyond current photonic quantum computing experimentation, thereby significantly reducing the likelihood of achieving near-term quantum supremacy in the context of Boson Sampling.

Equations

$$\mathcal{O}(n{2}^n+\text{poly}(m,n)),$$

$$\mathrm{Per}M=\frac{1}{2^{n-1}}\sum _{\delta }\left(\prod _{k=1}^n{\delta }_k\right)\prod _{j=1}^n\sum _{i=1}^n{\delta }_i{m}_{i,j}$$

$$q(\mathbf{z})=\frac{1}{\mu (\mathbf{z})}{\left|\mathrm{Per}{A}_{\mathbf{z}}\right|}^2\stackrel{\text{defn}}{=}\frac{1}{\mu (\mathbf{z})}{\left|\sum _{\sigma }\prod _{k=1}^n{a}_{z_k{\sigma }_k}\right|}^2,\mathbf{z}\in {\Phi }_{m,n},$$

$$p(\mathbf{r})=\frac{1}{n!}{\left|\mathrm{Per}{A}_{\mathbf{r}}\right|}^2=\frac{1}{n!}{\left|\sum _{\sigma }\prod _{i=1}^n{a}_{r_i{\sigma }_i}\right|}^2,\mathbf{r}\in [m{]}^n.$$

$$q(\mathbf{z})=\frac{n!}{\mu (\mathbf{z})}p(\mathbf{z})=\frac{n!}{\mu (\mathbf{z})}\frac{1}{n!}{\left|\mathrm{Per}{A}_{\mathbf{z}}\right|}^2=\frac{1}{\mu (\mathbf{z})}{\left|\mathrm{Per}{A}_{\mathbf{z}}\right|}^2,\mathbf{z}\in {\Phi }_{m,n},$$

$$p(r_1,\dots ,r_k)=\frac{(n-k)!}{n!}\sum _{c\in {\mathcal{C}}_k}{\left|\mathrm{Per}{A}_{r_1,\dots ,r_k}^c\right|}^2,k=1,\dots ,n,$$

$$n!\sum _{\mathbf{r}}p(\mathbf{r})$$

$$=\sum _{\mathbf{r}}\left(\sum _{\sigma }\prod _{i=1}^n{a}_{r_i{\sigma }_i}\right)\left(\sum _{\tau }\prod _{i=1}^n{\bar{a}}_{r_i{\tau }_i}\right),\tau \in \pi [n]$$

$$=\sum _{\sigma }\sum _{\tau }\left(\sum _{\mathbf{r}}\prod _{i=1}^n{a}_{r_i{\sigma }_i}{\bar{a}}_{r_i{\tau }_i}\right)$$

$$=\sum _{\sigma }\sum _{\tau }\prod _{i=1}^n\sum _{k=1}^m{a}_{k,{\sigma }_i}{\bar{a}}_{k{\tau }_i}.$$

$$n!p(r_1)$$

$$=\sum _{\sigma }\sum _{\tau }{a}_{r_1{\sigma }_1}{\bar{a}}_{r_1{\tau }_1}\prod _{i=2}^n\sum _{k=1}^m{a}_{k,{\sigma }_i}{\bar{a}}_{k{\tau }_i}.$$

$$n!p(r_1)=\sum _{\sigma }|a_{r_1{\sigma }_1}{|}^2=(n-1)!\sum _{k=1}^n|a_{r_1,k}{|}^2,$$

$$p(r_1)=\frac{1}{n}\left[|a_{r_1,1}{|}^2+\dots ,+|a_{r_1,n}{|}^2\right],r_1\in [m],$$

$$n!p(r_1,r_2)=\sum _{\sigma }\sum _{\tau }{a}_{r_1{\sigma }_1}{\bar{a}}_{r_1{\tau }_1}{a}_{r_2{\sigma }_2}{\bar{a}}_{r_2{\tau }_2}\left(\prod _{i=3}^n\sum _{k=1}^m{a}_{k,{\sigma }_i}{\bar{a}}_{k{\tau }_i}\right).$$

$$n!p(r_1,r_2)$$

$$=(n-2)!\sum _{c\in {\mathcal{C}}_2}{\left(\sum _{\sigma \in \pi (c)}{a}_{r_1{\sigma }_1}{\bar{a}}_{r_2{\sigma }_2}\right)}^2$$

$$=(n-2)!\sum _{c\in {\mathcal{C}}_2}{\left|\mathrm{Per}{A}_{r_1,r_2}^c\right|}^2,$$

$$n!p(r_1,\dots ,r_k)=(n-k)!\sum _{c\in {\mathcal{C}}_k}{\left|\mathrm{Per}{A}_{r_1,\dots ,r_k}^c\right|}^2,$$

$$p(r_1,\dots ,r_k)=\frac{(n-k)!}{n!}\sum _{c\in {\mathcal{C}}_k}{\left|\mathrm{Per}{A}_{r_1,\dots ,r_k}^c\right|}^2,k=1,\dots ,n.$$

$$p(\mathbf{r})=p(r_1)p(r_2|r_1)p(r_3|r_1,r_2)\dots p(r_n|r_1,r_2,\dots ,r_{n-1}).$$

$$\sum _{k=1}^{n}k{2}^k\left(\genfrac{}{}{0ex}{}{n}{k}\right)=\frac{2}{3}n{3}^n,$$

$$\mathrm{Per}B=\sum _{\ell =1}^k{b}_{k,\ell }\mathrm{Per}{B}_{k,\ell }^{\diamond },$$

$$\mathrm{Per}{B}_{k,\ell }^{\diamond }=\frac{1}{2^{k-2}}\sum _{\delta }\left(\prod _{i=1}^{k-1}{\delta }_i\right)\prod _{j\in [k]\setminus \ell }{v}_j(\delta ),\ell \in [k],$$

$$\varphi (r_1,\dots ,r_k|\mathbf{\alpha })=\frac{1}{k!}{\left|\mathrm{Per}{A}_{r_1,\dots ,r_k}^{[n]\setminus \{{\alpha }_{k+1},\dots ,{\alpha }_n\}}\right|}^2,k=1,\dots ,n-1.$$

$$p(r_1,\dots ,r_{n-1})=\frac{1}{n}\sum _{j=1}^n\frac{1}{(n-1)!}{\left|\mathrm{Per}{A}_{r_1,\dots ,r_{n-1}}^{[n]\setminus j}\right|}^2={\mathbb{E}}_{{\alpha }_n}\{\varphi (r_1,\dots ,r_{n-1}|\mathbf{\alpha })\},$$

$$\sum _{r_k}\varphi (r_1,\dots ,r_k|\mathbf{\alpha })={\mathbb{E}}_{{\alpha }_k|{\alpha }_{k+1},\dots ,{\alpha }_n}\{\varphi (r_1,\dots ,r_{k-1}|\mathbf{\alpha })\},k=2,\dots ,n-1,$$

$$\varphi (\mathbf{r}|\mathbf{\alpha })=\prod _{k=1}^n{e}_k/d_k=\varphi (r_1|\mathbf{\alpha })\varphi (r_2|r_1,\mathbf{\alpha })\varphi (r_3|r_1,r_2,\mathbf{\alpha })\dots \varphi (r_n|r_1,r_2,\dots ,r_{n-1},\mathbf{\alpha }),$$

$$\varphi (r_1|\mathbf{\alpha })=|\mathrm{Per}{A}_{r_1}^{[n]\setminus \{{\alpha }_2,\dots ,{\alpha }_n\}}{|}^2=|\mathrm{Per}{A}_{r_1}^{{\alpha }_1}{|}^2=|a_{r_1,{\alpha }_1}{|}^2,r_1\in [m].$$

$${\left|\mathrm{Per}{A}_{r_1,\dots ,r_k}^{{\alpha }_1,\dots ,{\alpha }_k}\right|}^2={\left|\sum _{\ell =1}^k{a}_{r_k,{\alpha }_{\ell }}\mathrm{Per}{B}_{k,\ell }^{\diamond }\right|}^2,r_k\in [m],$$

$$\mathrm{Per}B=\frac{1}{2^{k-1}}\sum _{\delta }\left(\prod _{i=1}^k{\delta }_i\right)\prod _{j\in [k]}{v}_j^{\ast }(\delta ),$$

$$\mathbb{E}(s_i)$$

$$\mathbb{E}(s_i(s_i-1))$$

$$P_D⩽\sum _{i=1}^m\mathbb{P}(s_i⩾2)⩽\frac{1}{2}\sum _{i=1}^m\mathbb{E}{s}_i(s_i-1)=2\sum _{i=1}^m\sum _{k<\ell }|a_{i,k}{a}_{i,\ell }{|}^2.$$

Full text

The fastest known exact classical algorithm for the standard Boson Sampling problem requires $\mathcal{O}((^{m+n-1}_{\hphantom{m+}n})\,n2^{n})$ time to produce samples for a system with input size $n$ and $m$ output modes, making it infeasible for anything but the smallest values of $n$ and $m$. We give an algorithm that is much faster, running in $\mathcal{O}(n2^{n}+\poly(m,n))$ time and $\mathcal{O}(m)$ additional space. The algorithm is simple to implement and has low constant factor overheads. As a consequence our classical algorithm is able to solve the exact Boson Sampling problem for system sizes far beyond current photonic quantum computing experimentation, thereby significantly reducing the likelihood of achieving near-term quantum supremacy in the context of Boson Sampling.

1 Introduction

The promise of significantly faster quantum algorithms for problems of practical interest is one of the most exciting prospects in computer science. Most famously, a quantum computer would allow us to factorise integers in polynomial time (Shor,, 1997) and perform unstructured search on an input of $n$ elements in $\mathcal{O}(\sqrt{n})$ time (Grover,, 1996). In the short term however, despite impressive progress in recent years, there are still considerable challenges to overcome before we can build a fully universal quantum computer. Until this time, the question of how to design quantum experiments which fall short of fully universal computation but which still show significant quantum speed-up over known classical (that is non-quantum) algorithms remains of great interest. In a breakthrough contribution in STOC ’11 Aaronson and Arkhipov, (2011) proposed an experimental set-up in linear optics known as Boson Sampling as a potentially practical way to do just this.

Since its introduction, Boson Sampling has attracted a great deal of attention with numerous experimental efforts around the world attempting implementations for various problem sizes (see e.g. Bentivegna et al.,, 2015; Spring et al.,, 2013; Broome et al.,, 2013; Tillmann et al.,, 2013; Crespi et al.,, 2013; Spagnolo et al.,, 2014; Latmiral et al.,, 2016). The ultimate goal is to exhibit a physical quantum experiment of such a scale that it would be hard if not impossible to simulate the output classically and thereby to establish so-called ‘quantum supremacy’. In terms of the physical Boson Sampling experiment, $n$ corresponds to the number of photons and $m$ the number of output modes and increasing either of these is difficult in practice. Progress has therefore been slow (see Lund et al.,, 2017, and the references therein) with the current experimental record being $n=5,m=9$ (Wang et al.,, 2016).

Translated into conventional computing terms, the task is to generate independent random samples from a specific probability distribution on multisets of size $n$ with elements in the range $1$ through $m$. For Boson Sampling, the probability of each multiset is related to the permanent of an $n\times n$ matrix built from possibly repeated rows of a larger $m\times n$ matrix where the multiset determines the rows used in the construction. The Boson Sampling problem is to produce such samples either with a quantum photonic device or with classical computing hardware.

On the face of it, Boson Sampling is a challenging problem for classical computation. Aaronson and Arkhipov, (2011) showed that exact Boson Sampling is not solvable in polynomial time by a classical computer unless $\P^{\#\P}=\BPP^{\NP}$ and the polynomial hierarchy collapses to the third level.

Without a clear understanding of the classical complexity of the problem it is difficult to determine the minimum experimental size needed to establish quantum supremacy. In the absence of a fast classical algorithm for Boson Sampling, speculation about the minimum size of $n$ (with $m=n^{2}$) has varied a great deal in the literature from e.g. $7$ (Latmiral et al.,, 2016), between $20$ and $30$ (Aaronson and Arkhipov,, 2014) and all the way up to $50$ (Lund et al.,, 2017). Broadly speaking the lowest estimate corresponds to the limits of the previous fastest classical Boson Sampling algorithm and the highest with an assumption that the polynomial hierarchy will not collapse.

Brute force evaluation of the probabilities of each multiset as a preliminary to random sampling, requires the calculation of $(^{m+n-1}_{\hphantom{m+}n})$ permanents of $n\times n$ matrices, each one of which takes $\mathcal{O}(n2^{n})$ time with the fastest known algorithms. If $m\geqslant n^{2}$, as considered by Aaronson and Arkhipov, (2011), the total running time is $\Theta(2^{n}\mathrm{e}^{n}(m/n)^{n}n^{1/2})$. The space usage of such a naive approach can be reduced to the storage needed for the required sample size (see e.g. Efraimidis,, 2015), however the running time for anything but the smallest values of $n$ and $m$ remains prohibitive. Previous to the work we present here, no faster provably correct classical Boson Sampling algorithm was known.

Recently Neville et al., (2017) gave strong numerical evidence suggesting that a specially designed Markov chain Monte Carlo (MCMC) implementation can sample approximately from the collision-free version of Boson Sampling. The time cost per sample for their MCMC method for problems with size $n\leqslant 30$ corresponds to that of computing around two hundred $n\times n$ permanents, making it computationally feasible for the problem sizes they consider. This was the first practical evidence that classical Boson Sampling might be possible faster than the naive brute force approach and is a significant motivation for the work we present here.

We show that Boson Sampling on a classical computer can be performed both exactly and dramatically faster than the naive algorithm, significantly raising the bar for the minimum size of quantum experiment needed to establish quantum supremacy. Our sampling algorithm also provides the probability associated with any given sample at no extra computational cost.

Theorem 1**.**

The time complexity of Boson Sampling is:

[TABLE]

where $\poly(m,n)=\mathcal{O}(mn^{2}).$ The additional space complexity on top of that needed to store the input is $\mathcal{O}(m)$.

A particularly attractive property of the running time we give is that the exponentially growing term depends only on $n$, the size of the multisets we wish to sample, and not $m$ which is potentially much larger. In the most natural parameter regime when $m$ is polynomial in $n$ our algorithm is therefore approximately $(^{m+n-1}_{\hphantom{m+}n})$ times faster than the fastest previous method. Our algorithm is both straightforward to implement and has small constant factor overheads, with the total computational cost to take one sample approximately equivalent to that of calculating two $n\times n$ permanents. This gives a remarkable resolution to the previously open question of the relationship between the classical complexity of computing the permanent of an $n\times n$ matrix and Boson Sampling. For an implementation of the algorithm, see Clifford and Clifford, (2017).

Based on the time complexity we demonstrate in Theorem 1 and bearing in mind recent feasibility studies for the calculation of large permanents (Wu et al.,, 2016) we find support for $n=50$ as the threshold for quantum supremacy for the exact Boson Sampling problem.

Our sampling algorithm also applies directly to the closely related problem of scattershot Boson Sampling (Bentivegna et al.,, 2015). From a mathematical perspective this produces no additional complications beyond the specification of a different $m\times n$ matrix. Once the matrix is specified we can apply our new sampling algorithm directly.

Our techniques draw heavily on the highly developed methods of hierarchical sampling commonly used in Bayesian computation (see e.g. Liu,, 2008; Green et al.,, 2015); specially the use of auxiliary variables and conditional simulation. We believe this cross-fertilisation of ideas is novel and opens up the development of classical algorithms for sampling problems in quantum computation. It raises the intriguing question of whether similar techniques can be applied to get significant improvements in time and space for other proposed quantum supremacy projects such as, for example, the problem of classically simulating quantum circuits (Boixo et al.,, 2016).

The algorithmic speed-up we achieve relies on a number of key innovations. As a first step we expand the sample space of multisets into a product set $[m]^{n}$, i.e. the space of all arrays $\mathbf{r}=(r_{1},\dots,r_{n})$, with each element in the range $1$ to $m$. This simple transformation allows us to expose combinatorial properties of the sampling problem which were otherwise opaque. We then develop explicit expressions for the probabilities of subsequences of arrays in this set. Application of the chain rule for conditional probabilities then leads to our first speed-up with an $\mathcal{O}(mn3^{n})$ time algorithm for Boson Sampling. For our main result in Theorem 1, we introduce two further innovations. First we exploit the Laplace expansion of a permanent to give an efficient amortised algorithm for the computation of the permanent of a number of closely related matrices. We then expand the sample space further with an auxiliary array representing column indices and show how a conditioning argument allows Boson Sampling to be viewed as sampling from a succession of coupled conditional distributions. The final result is our improved time complexity of $\mathcal{O}(n2^{n}+mn^{2})$.

A consequence of our analysis is that we are also able to provide an efficiently computable bound on the probability that a sampled multiset drawn from the Boson Sampling distribution will contain duplicated elements. This enables us to bound the performance of rejection sampling variants of our sampling algorithm when applied to the so-called ‘collision-free’ setting described by Aaronson and Arkhipov,.

One further application of our new sampling algorithm lies in tackling the question of whether the output of physical Boson Sampling devices correspond to theoretical predictions. There has been considerable interest in this area (Aaronson and Arkhipov,, 2014; Tichy et al.,, 2014; Wang and Duan,, 2016; Walschaers et al.,, 2016) however previously developed methods have had to work under the assumption that we cannot sample from the Boson Sampling distribution for even moderate sizes of $n$ and $m$. Armed with our new sampling algorithm, we can now sample directly from the Boson Sampling distribution for a much larger range of parameters, allowing a far greater range of statistical tests to be applied.

Related work and computation of the permanent

Terhal and DiVincenzo, (2004) were the first to recognise that studying the complexity of sampling from low-depth quantum circuits could be useful for demonstrating a separation between quantum and classical computation. Since that time the goal of finding complexity separations for quantum sampling problems has received a great deal of interest. We refer the interested reader to Lund et al., (2017) and Harrow and Montanaro, (2017) for recent surveys on the topic.

There is an extensive literature on the problems of sampling from unwieldy distributions with applications in all the main scientific disciplines. Much of this work has relied on MCMC methods where the target distribution can be identified as the equilibrium distribution of a Markov chain (see e.g. Hastings,, 1970). For Boson Sampling the most closely related work to ours is that of Neville et al., (2017) mentioned earlier. In their paper the authors also set out the hurdles that need to be overcome for a physical Boson Sampling experiment to reach higher values of $n$. Their MCMC approach does not however permit them to estimate how well their sampling method approximates the true Boson Sampling distribution, other than through numerical experimentation. Our new exact sampling algorithm is therefore the first provably correct algorithmic speed-up for Boson Sampling. In fact our algorithm also has smaller constant factor overheads, costing approximately two permanent calculations of $n\times n$ matrices for all values of $n$ (assuming $m$ remains polynomial in $n$).

There are MCMC algorithms which permit exact sampling, for example via ‘coupling from the past’, where guarantees of performance can be provided (Propp and Wilson,, 1996). However, the range of problems for which this approach has proven applicable is necessarily more limited and there has yet to be a successful demonstration for the Boson Sampling problem.

Permanents:

Computation of the permanent of a matrix was shown to be #¶-hard by Valiant, (1979) and hence is unlikely to have a polynomial time solution. Previously Ryser, (1963) had proposed an $O(n^{2}2^{n})$ time algorithm to compute the permanent of an $n\times n$ matrix. This was sped up to $O(n2^{n})$ time by Nijenhuis and Wilf, (1978) and many years later a related formula with the same time complexity was given by Glynn, (2010). It is Glynn’s computational formula for the permanent which we will use as the basis of our sampling algorithms.

Let $M=(m_{i,j})$ be an $n\times n$ matrix with $m_{i,j}\in\mathbb{C}$, and let $\mathop{\textup{Per}}M$ denote its permanent, then

[TABLE]

where $\delta\in\{-1,1\}^{n}$ with $\delta_{1}=1$.

A naive implementation of this formula would require $\mathcal{O}(n^{2}2^{n})$ time however by considering the $\delta$ arrays in Gray code order this is reduced to $O(n2^{n})$ time. Note that further reduction in the time complexity may be possible when $M$ has repeated rows or columns; see Tichy, (2011, Appendix B) and Shchesnovich, (2013, Appendix D). In the special case of computation of the permanent over $\poly(n)$ bit integers, Björklund, (2016) has a faster $\mathcal{O}(n2^{n-\Omega(\sqrt{n/\log{n}})})$ time algorithm for computing the permanent.

2 Problem specification

In mathematical terms the Boson Sampling problem is as follows.

Let $m$ and $n$ be positive integers and consider all multisets of size $n$ with elements in $[m]$, where $[m]=\{1,\dots,m\}$. Each multiset can be represented by an array $\mathbf{z}=(z_{1},\dots,z_{n})$ consisting of elements of the multiset in non-decreasing order: $z_{1}\leqslant z_{2}\leqslant\dots\leqslant z_{n}$. We denote the set of distinct values of $\mathbf{z}$ by $\Phi_{m,n}$. The cardinality of $\Phi_{m,n}$ is known to be $(^{m+n-1}_{\hphantom{m+}n})$ – see Feller, (1968), for example. We define $\mu(\mathbf{z})=\prod_{j=1}^{m}s_{j}!$ where $s_{j}$ is the multiplicity of the value $j$ in $\mathbf{z}$.

Now let $A=(a_{ij})$ be a complex valued $m\times n$ matrix consisting of the first $n$ columns of a given $m\times m$ Haar random unitary matrix. For each $\mathbf{z}$, build an $n\times n$ matrix $A_{\mathbf{z}}$ where the $k$-th row of $A_{\mathbf{z}}$ is row $z_{k}$ in $A$ for $k=1,\dots,n$ and define a probability mass function (pmf) on $\Phi_{m,n}$ as

[TABLE]

where $\mathop{\textup{Per}}A_{\mathbf{z}}$ is the permanent of $A_{\mathbf{z}}$ and in the definition the summation is for all $\sigma\in\pi[n]$, the set of permutations of $[n]$.

The Cauchy-Binet formula for permanents (see Marcus and Minc,, 1965, page 579) can be used to show the function $q$ is indeed a pmf. There is a demonstration from first principles in the proof of Lemma 1. See also Aaronson and Arkhipov, (2011, Theorem 3.10) for a demonstration from physical principles.

The computing task is to simulate random samples from the pmf $q(\mathbf{z})$.

3 Exact Boson Sampling Algorithms

We give two algorithms for Boson Sampling. Algorithm A runs in $\mathcal{O}(mn3^{n})$ time to create a single sample. This is already a significant speed-up over the $\Theta{\left((^{m+n-1}_{\hphantom{m+}n})n2^{n}\right)}$ time complexity of the fastest solution in the literature. Algorithm A prepares the way for Algorithm B our main result with running time specified in Theorem 1.

We will repeatedly need to sample from an unnormalised discrete probability mass function (weight function). To establish the complexity of our Boson Sampling algorithms we need only employ the most naive linear time method although faster methods do exist. In particular, if the probability mass function is described by an array of length $m$ then it is possible to sample in $\mathcal{O}(1)$ time after $\mathcal{O}(m)$ preprocessing time, using Walker’s alias method (Walker,, 1974; Kronmal and Peterson Jr,, 1979).

Expanding the sample space:

We start by expanding the sample space $\Phi_{m,n}$ and consider a distribution on the product space $[m]^{n}$, i.e. the space of all arrays $\mathbf{r}=(r_{1},\dots,r_{n})$, with each element in the range $1$ to $m$. Since the permanent of a matrix is invariant to row reordering, we have $\mathop{\textup{Per}}A_{\mathbf{r}}=\mathop{\textup{Per}}A_{\mathbf{z}}$ when $\mathbf{r}$ is any permutation of $\mathbf{z}$. Furthermore the number of distinguishable rearrangements of $\mathbf{z}$ is $n!/\mu(\mathbf{z})$.

We claim that in order to sample from $q(\mathbf{z})$ we can equivalently sample from the pmf

[TABLE]

In other words sample $\mathbf{r}$ from $p(\mathbf{r})$ and then rearrange the elements of $\mathbf{r}$ in non-decreasing order to give $\mathbf{z}$.

This follows since for any $\mathbf{z}$ there are $n!/\mu(\mathbf{z})$ equally likely values of $\mathbf{r}$ in the expanded sample space, i.e. $p(\mathbf{r})=p(\mathbf{z})$ for all such values of $\mathbf{r}$. So by the addition rule for probabilities

[TABLE]

as claimed.

This straightforward reformulation will make the task of efficient Boson Sampling significantly simpler. Our approach focuses on the pmf $p(\mathbf{r})$.

We begin by deriving the marginal pmfs of the leading subsequences of $(r_{1},\dots,r_{n})$ in the following lemma. This will then provide the pmfs of successive elements of $(r_{1},\dots,r_{n})$ conditional on previous values. In that way we are able to rewrite $p(r_{1},\dots,r_{n})$ as a product of conditional pmfs, using the chain rule of probability, and exploit the chain for progressively simulating $r_{1},\dots,r_{n}$.

Lemma 1**.**

The joint pmf of the subsequence $(r_{1},\dots,r_{k})$ is given by

[TABLE]

where $\mathcal{C}_{k}$ is the set of $k$-combinations taken without replacement from $[n]$ and $A_{r_{1},\dots,r_{k}}^{c}$ is the matrix formed from rows $(r_{1},\dots,r_{k})$ of the columns $c$ of $A$.

Proof.

The convention will be that $c$ is a set and $(r_{1},\dots,r_{k})$ is an array with possibly multiple instances of its elements. The invariance of permanents under column and row permutations removes any ambiguity in the interpretation of $\mathop{\textup{Per}}A^{c}_{r_{1},\dots,r_{k}}$.

For completeness we start by confirming that $p(\mathbf{r})$ is a pmf, in other words $\sum_{\mathbf{r}}p(\mathbf{r})=1.$ From the definition of $p(\mathbf{r})$ and multiplying by $n!$, we have

[TABLE]

The product is [math] when $\sigma_{i}\neq\tau_{i}$ for any $i$, by the orthonormal property of $A$. Otherwise when $\sigma=\tau$ the product is $1$ so the final expression reduces to $n!$ and $\sum_{\mathbf{r}}p(\mathbf{r})=1$ as expected.

Now looking at the case $k=1$ we have similarly

[TABLE]

Again for this to be non-zero, each of the terms in the product must be $1$. This means $\sigma_{i}=\tau_{i},i=2,\dots,n$ which implies $\sigma=\tau$.

It follows that

[TABLE]

as claimed, so that

[TABLE]

the sum of the squared moduli of elements of row $r_{1}$ divided by $n$.

In the same way,

[TABLE]

Again the only non-zero case is when $\sigma_{i}=\tau_{i},i=3,\dots,n$. However it does not imply that $\sigma=\tau$, since we can also have $\sigma_{1}=\tau_{2},\sigma_{2}=\tau_{1}$.

Consequently,

[TABLE]

where $\mathcal{C}_{2}$ is the set of $2$-combinations (without replacement) from $[n]$, $\pi(c)$ is the set of permutations of $c$ and $A_{r_{1},r_{2}}^{c}$ is a matrix formed from rows $(r_{1},r_{2})$ and columns $c$ of $A$.

In the same way

[TABLE]

so that

[TABLE]

∎

Our first Boson Sampling algorithm which we term Algorithm A samples from $p(\mathbf{r})$ by expressing this as a chain of conditional pmfs,

[TABLE]

Sample $r_{1}$ from $p(r_{1})$, $r_{1}\in[m]$. Subsequently for stages $k=2,\dots,n$, sample $r_{k}$ from the conditional pmf $p(r_{k}|r_{1},\dots,r_{k-1})$, $r_{k}\in[m]$ with $(r_{1},\dots,r_{k-1})$ fixed. At completion, the array $(r_{1},\dots,r_{n})$ at stage $n$ will have been sampled from the pmf $p(\mathbf{r})$. Sorting $(r_{1},\dots,r_{n})$ in non-decreasing order, the resulting multiset $\mathbf{z}$ is sampled from the Boson Sampling distribution $q(\mathbf{z})$.

Theorem 2**.**

Algorithm A samples from the Boson Sampling distribution with time complexity $\mathcal{O}(mn3^{n})$ per sample. The additional space complexity of Algorithm A on top of that needed to store the input is $\mathcal{O}(m)$ words.

Proof.

The correctness follows from the chain rule for conditional probabilities and so we now analyse the running time. From Lemma 1, evaluation of the pmf $p(r_{1})$ for $r_{1}\in[m]$ involves the sum of squared moduli in each row of $A$, a total of $\mathcal{O}(mn)$ operations. Sampling from this pmf takes $\mathcal{O}(m)$ time. More generally we can write $p(r_{k}|r_{1},\dots,r_{k-1})=p(r_{1},\dots,r_{k})/p(r_{1},\dots,r_{k-1})$. Note that $r_{k}$ does not appear in the denominator, so to sample $r_{k}$ from this conditional pmf we can equivalently sample from the pmf that is proportional to the numerator, considered purely as function of $r_{k}$ with $r_{1},\dots,r_{k-1}$ fixed. The numerator involves the calculating of several $k\times k$ permanents. With the best known algorithms for computing the permanent, each requires $\mathcal{O}(k2^{k})$ operations. There are $(^{n}_{k})$ terms in $\mathcal{C}_{k}$ for each value of $k$, making a combined operation count at stage $k$ of $\mathcal{O}(m(^{n}_{k})k2^{k})$. Using

[TABLE]

the total operation count is then $\mathcal{O}(mn3^{n})$. We only need to store the result of one permanent calculation at a time and also one list of $m$ probabilities at a time. This gives a total space usage of $\mathcal{O}(m)$ excluding the space needed to store the original matrix $A$. ∎

We now prepare to prove Theorem 1.

Laplace expansion:

We make extensive use of the Laplace expansion for permanents (see Marcus and Minc,, 1965, page 578), namely that for any $k\times k$ matrix $B=(b_{i,j})$,

[TABLE]

where $B^{\diamond}_{k,\ell}$ is the submatrix with row $k$ and column $\ell$ removed. Note that $B^{\diamond}_{k,\ell}$ only depends on $B^{\diamond}_{k}$ the submatrix of $B$ with the $k$-th row removed. An important consequence is that when $B$ is modified by changing its $k$-th row, the modified permanent can be calculated in $\mathcal{O}(k)$ steps, provided the values $\{\mathop{\textup{Per}}B^{\diamond}_{k,\ell}\}$ are available.

We show that computation of all of the values $\{\mathop{\textup{Per}}B^{\diamond}_{k,\ell},\ell\in[k]\}$ has the same time complexity as computing $\mathop{\textup{Per}}B$, the permanent of a single $k\times k$ matrix. We will later take advantage of this fast amortised algorithm in the proof of Theorem 1 to quickly compute a set of permanents of matrices each with one row differing from the other.

Lemma 2**.**

Let $B$ be a $k\times k$ complex matrix and let $\{B^{\diamond}_{k,\ell}\}$ be submatrices of $B$ with row $k$ and column $\ell$ removed, $\ell\in[k]$. The collection $\{\mathop{\textup{Per}}B^{\diamond}_{k,\ell},\ell\in[k]\}$ can be evaluated jointly in $\mathcal{O}(k2^{k})$ time and $\mathcal{O}(k)$ additional space.

Proof of Lemma.

By applying Glynn’s formula (1) for a given value of $\ell$ we have:

[TABLE]

where $\delta\in\{-1,1\}^{k-1}$ with $\delta_{1}=1$ and $v_{j}(\delta)=\sum_{i=1}^{k-1}\delta_{i}b_{ij}$.

Exploiting the usual trick of working through values of $\delta$ in Gray code order, the terms $\{v_{j}(\delta),j\in[k]\}$ can be evaluated in $\mathcal{O}(k)$ combined time for every new $\delta$. This is because successive $\delta$ arrays will differ in one element. The product of the $v_{j}(\delta)$ terms can therefore also be computed in $\mathcal{O}(k)$ time giving $\mathcal{O}(k2^{k})$ time to compute $\mathop{\textup{Per}}B^{\diamond}_{k,\ell}$ for a single value of $\ell$, but of course this has to be replicated $k$ times to cover all values of $\ell$.

To compute $\{\mathop{\textup{Per}}B^{\diamond}_{k,\ell}\ell\in[k]\}$ more efficiently we observe that each product $\prod_{j\in[k]\setminus\ell}v_{j}(\delta)$ can be expressed as $f_{\ell}b_{\ell}$ where $f_{\ell}=\prod_{j=1}^{\ell-1}v_{j}(\delta),\ell=2,\dots,k$ and $b_{\ell}=\prod_{j=\ell+1}^{k}v_{j}(\delta),\ell=1,\dots,k-1$ are forward and backward cumulative products, with $f_{1}=b_{k}=1$. We can therefore compute all of the partial products $\prod_{j\in[k]\setminus\ell}v_{j}(\delta)$ in $\mathcal{O}(k)$ time, giving an overall total time complexity of $\mathcal{O}(k2^{k})$ for jointly computing $\{\mathop{\textup{Per}}B^{\diamond}_{k,\ell},\ell\in[k]\}$. Other than the original matrix, space used is dominated by the array of $k$ accumulating partial products. ∎

Furthermore the computation time has constant factor overheads similar to that of computing $\mathop{\textup{Per}}B$ (see Appendix).

Auxiliary column indices:

We now expand the sample space even further with an auxiliary array $\boldsymbol{\alpha}=(\alpha_{1},\dots,\alpha_{n})$, where $\boldsymbol{\alpha}\in\pi[n]$, the set of permutations of $[n]$.

As with Lemma 1, the purpose is to create a succession of pmfs for leading subsequences of the array $(r_{1},\dots,r_{n})$ and thereby provide successive conditional pmfs so that $r_{1},\dots,r_{n}$ can be simulated progressively. The novelty here is to introduce an overall conditioning variable $\boldsymbol{\alpha}$ so that sampling can be carried out more rapidly under that condition while demonstrating that the marginal distribution still correctly provides the pmf $p(\mathbf{r})$ in (3).

Define

[TABLE]

For notational convenience let $e_{k}=\phi(r_{1},\dots,r_{k}|\boldsymbol{\alpha})$ and $d_{k}=\sum_{r_{k}}\phi(r_{1},\dots,r_{k}|\boldsymbol{\alpha})$ for $k=1,\dots,n-1$ with $e_{n}=p(r_{1},\dots,r_{n})$ and $d_{n}=p(r_{1},\dots,r_{n-1})$. Note that $\phi(r_{1},\dots,r_{k}|\boldsymbol{\alpha})$ is a pmf on $[m]^{k}$, equivalent to (3) using the subset $\{\alpha_{1},\dots,\alpha_{k}\}$ of columns of $A$, so that in particular $d_{1}=1$.

Lemma 3**.**

With the preceding notation, let $\phi(\mathbf{r}|\boldsymbol{\alpha})=\prod_{k=1}^{n}e_{k}/d_{k}$ then $p(\mathbf{r})=\mathbb{E}_{\boldsymbol{\alpha}}\{\phi(\mathbf{r}|\boldsymbol{\alpha})\}$ where the expectation is taken over $\boldsymbol{\alpha}$, uniformly distributed on $\pi[n]$ for fixed $\mathbf{r}$.

Proof.

For the particular case $k=n-1$, Lemma 1 shows that $p(r_{1},\dots,r_{n-1})$ can be written as a mixture of component pmfs and expressed as an expectation over $\alpha_{n}$, i.e.

[TABLE]

where $\alpha_{n}$ is uniformly distributed on $[n]$. In other words $d_{n}$ is the expectation of $e_{n-1}$ over $\alpha_{n}$. More generally when the leading subsequence is $(r_{1},\dots,r_{k})$ and the active columns of $A$ are $[n]\setminus\{\alpha_{k+1},\dots,\alpha_{n}\}$, we have

[TABLE]

where the expectation is taken over $\alpha_{k}$ uniformly distributed on $[n]\setminus\{\alpha_{k+1},\dots,\alpha_{n}\}$ with $(\alpha_{k+1},\dots,\alpha_{n})$ and $\mathbf{r}$ fixed. This means $d_{k}$ is the conditional expectation of $e_{k-1}$.

We now rewrite $\phi(\mathbf{r}|\boldsymbol{\alpha})$ as $e_{n}\prod_{k=2}^{n}e_{k-1}/d_{k}$ since $d_{1}=1$ and start by considering the expectations of the terms $e_{k-1}/d_{k},k=2,\dots,n$ with respect to $\alpha_{2}$ given $\alpha_{3},\dots,\alpha_{n}$. With the exception of $e_{1}$ these terms are specified by subsets of $\{\alpha_{3},\dots,\alpha_{n}\}$, so they remain fixed under the conditioning. Since the conditional expectation of $e_{1}$ is $d_{2}$, the conditional expectation of $e_{n}\prod_{k=2}^{n}e_{k-1}/d_{k}$ is $e_{n}\prod_{k=3}^{n}e_{k-1}/d_{k}$, i.e. there is one fewer term in the product. Proceeding in this way and successively taking conditional expectations with respect to $\alpha_{k-1}$ uniformly distributed on $[n]\setminus\{\alpha_{k},\dots,\alpha_{n}\}$ conditional on $(\alpha_{k},\dots,\alpha_{n})$, the product telescopes to $e_{n}(e_{n-1}/d_{n})=p(r_{1},\dots,r_{n})e_{n-1}/d_{n}$. Finally taking expectations over $\alpha_{n}$, this reduces to $p(r_{1},\dots,r_{n})$, since $d_{n}$ is the expected value of $e_{n-1}$. Note a random permutation $(\alpha_{1},\dots,\alpha_{n})$ is generated starting from $\alpha_{n}$ in reverse order by successively selecting elements for inclusion at random from what remains. By the chain rule of expectation, we have then shown that $\mathbb{E}_{\boldsymbol{\alpha}}\{\phi(\mathbf{r}|\boldsymbol{\alpha})\}=p(\mathbf{r})$ as claimed. ∎

We can now describe our second and faster Boson Sampling algorithm. This algorithm is the basis of our main result, Theorem 1. We start by sampling a random permutation $\boldsymbol{\alpha}=(\alpha_{1},\dots,\alpha_{n})$ uniformly distributed on $\pi[n]$. As in Algorithm A we consider a chain of conditional pmfs,

[TABLE]

where now $\phi(r_{k}|r_{1},\dots,r_{k-1},\boldsymbol{\alpha})=e_{k}/d_{k}$, with the notation defined prior to Lemma 3.

At stage 1, we sample $r_{1}$ from the pmf $\phi(r_{1}|\boldsymbol{\alpha})$. Subsequently for stages $k=2,\dots,n$, we sample $r_{k}$ from $\phi(r_{k}|r_{1},\dots,r_{k-1},\boldsymbol{\alpha})$ with $(r_{1},\dots,r_{k-1})$ fixed, exploiting the Laplace expansion for permanents. At completion, the array $(r_{1},\dots,r_{n})$ at stage $n$ will have been sampled from the pmf $p(\mathbf{r})$. Sorting $(r_{1},\dots,r_{n})$ in non-decreasing order, the resulting multiset $\mathbf{z}$ is sampled from the Boson Sampling distribution $q(\mathbf{z})$.

Proof of Theorem 1.

It follows from the chain rule for conditional probabilities that for a given fixed $\boldsymbol{\alpha}$ Algorithm B samples from $\phi(\mathbf{r}|\boldsymbol{\alpha})$. As $\boldsymbol{\alpha}$ is uniformly distributed on $\pi[n]$ the algorithm therefore samples from the marginal distribution $\mathbb{E}_{\boldsymbol{\alpha}}\{\phi(\mathbf{r}|\boldsymbol{\alpha})\}$ which by Lemma 3 is the pmf $p(\mathbf{r})$ in (3). The complexity is established as follows. At stage 1, the pmf to be sampled is

[TABLE]

Sampling takes $\mathcal{O}(m)$ operations. At stage 2 for fixed $r_{1}$, we need to sample $r_{2}$ from the pmf $\phi(r_{2}|r_{1},\boldsymbol{\alpha})$. But since $\phi(r_{2}|r_{1},\boldsymbol{\alpha})=e_{2}/d_{2}$ and $d_{2}$ does not involve $r_{2}$, we can sample from the pmf that is proportional to $e_{2}$ considered simply as a function of $r_{2}$ with $r_{1}$ fixed. Since $e_{2}$ is proportional to $|\mathop{\textup{Per}}A^{\alpha_{1},\alpha_{2}}_{r_{1},r_{2}}|^{2}$ this involves calculating $m$ permanents of $2\times 2$ matrices and then sampling; a further $\mathcal{O}(m)$ operations.

At stage $k$ with $r_{1},\dots,r_{k-1}$ already determined, we need to sample $r_{k}$ from the pmf proportional to $|\mathop{\textup{Per}}A^{\alpha_{1},\dots,\alpha_{k}}_{r_{1},\dots,r_{k}}|^{2}$ considered simply as a function of $r_{k}$. Let $B=A^{\alpha_{1},\dots,\alpha_{k}}_{r_{1},\dots,r_{k}}$ for some arbitrary value of $r_{k}$, i.e. $b_{i,j}=a_{r_{i},\alpha_{j}}$ for $i,j=1,\dots,k$. Using the Laplace expansion and taking squared moduli

[TABLE]

where $B^{\diamond}_{k,\ell}$ is defined in the statement of Lemma 2. From Lemma 2 all the values $\{\mathop{\textup{Per}}B^{\diamond}_{k,\ell}\}$ can be computed in $\mathcal{O}(k2^{k})$ time. An array of length $m$ is then obtained from (4) by summing $k$ terms for each $r_{k}\in[m]$. Taking a single sample from the pmf proportional to the array takes $\mathcal{O}(m)$ time. This gives a total time complexity for stage $k$ of $\mathcal{O}(k2^{k})$ + $\mathcal{O}(mk)$.

The total time for all stages $k=1,\dots,n$ is therefore $\mathcal{O}(n2^{n})+\mathcal{O}(mn^{2})$. In practical terms, the exponentially growing term is approximately equivalent to the cost of evaluating two $n\times n$ permanents. The space usage is dominated by storing a single unnormalised pmf of size $\mathcal{O}(m)$ prior to sampling. ∎

Remarks:

At the final stage when $r_{n}$ is selected the value of $\left|\mathop{\textup{Per}}A^{\alpha_{1},\dots,\alpha_{n}}_{r_{1},\dots,r_{n}}\right|^{2}=\left|\mathop{\textup{Per}}A_{\mathbf{r}}\right|^{2}$ will have been computed as in (4). Hence $q(\mathbf{r})$ the pmf of the sample $\mathbf{r}$ is available at no extra computational cost. Finally note that when only a single value is to be sampled from $q(\mathbf{r})$ and since the columns of a Haar random unitary matrix are already in a random order, the random permutation $(\alpha_{1},\dots,\alpha_{n})$ in Algorithm B can be taken to be $(1,2,\dots,n)$

4 Acknowledgements

RC was funded by EPSRC Fellowship EP/J019283/1. Thank you to Ashley Montanaro, Alex Neville, Anthony Laing and members of the Bristol QET Labs for introducing us to the Boson Sampling problem and for a number of detailed and helpful conversations.

Appendix

Constant factor overheads in Lemma 2:

To compare the constant factor overheads in Lemma 2, we start by looking at $\mathop{\textup{Per}}B$. From Glynn’s formula we have

[TABLE]

where $\delta\in\{-1,1\}^{k}$ with $\delta_{1}=1$ and $v^{*}_{j}(\delta)=\sum_{i=1}^{k}\delta_{i}b_{ij}$. As we work through the values of $\delta$ in Gray code order, updating $v^{*}_{j}(\delta),j\in[k]$ requires $k$ additions. Evaluating the product requires $k-1$ multiplications and there is one further addition as the sum is accumulated. There are $2^{k-1}$ steps making a total of $2k2^{k-1}=k2^{k}$ operations.

For the calculation of $\{\mathop{\textup{Per}}B^{\diamond}_{k,\ell},\ell\in[k]\}$ in Lemma 2, at each stage there are $k$ additions for updating $\{v_{j}(\delta),j\in[k]\}$, then $k-2$ multiplications for calculating each of the sets $\{f_{\ell}\}$, $\{b_{\ell}\}$, and $\{f_{\ell}b_{\ell}\}$. Having thereby obtained all the partial products $\prod_{j\in[k]}v_{j}(\delta)/v_{\ell}(\delta)$, a further $k$ additions are required as the $k$ partial sums are accumulated. So each of the $2^{k-2}$ steps requires $5k-6$ operations. Consequently the operation count is approximately 25% larger than that of $\mathop{\textup{Per}}B$. In practice addition and multiplication with complex arithmetic will have different costs in terms of floating point operations.

Collision probability:

Suppose that $\mathbf{z}$ is sampled from $q(\mathbf{z})$ in (2) and that $\mathbf{s}$ is the associated array of counts. The probability of duplicate elements in $\mathbf{z}$ is the probability that one of the elements of $\mathbf{s}$ is larger than $1$. In general this probability will depend on the underlying Haar random unitary $A$.

From Lemma 1 we have

[TABLE]

Let us denote the $A$-specific duplication probability, i.e. the chance of seeing duplicated values in $\mathbf{z}$, by $P_{D}$. We have from (Collision probability:)

[TABLE]

Note that the expected value of the upper bound is $n(n-1)/(m+1)\sim n^{2}/m$ when averaged over $A$ from the Haar random unitary class. See Aaronson and Arkhipov, (2011, Theorem C.4) and Arkhipov and Kuperberg, (2012) for comparison.