Simulating Gaussian boson sampling on graphs in polynomial time
November 20, 2025
Abstract
We show that a distribution related to Gaussian Boson Sampling (GBS) on graphs can be sampled classically in polynomial time. Graphical applications of GBS typically sample from this distribution, and thus quantum algorithms do not provide exponential speedup for these applications. We also show that another distribution related to Boson sampling can be sampled classically in polynomial time.
1 Introduction↩︎
Bosons are subatomic particles with quantum spins. Boson sampling (BS) is a model of quantum computation based on interacting bosons that can be implemented using an optical network. It was proposed by Aaronson and Arkhipov [1] as a model that may have quantum advantage in the sense that sampling approximately from the output distribution can be done on a Boson computer. However, assuming two conjectures - the permanent anti-concentration conjecture (concerning the permanent of a matrix of i.i.d.Gaussians) - and the permanent of Gaussian conjecture (that estimating this permanent is #P-hard), they show that classically sampling from the same distribution would collapse the polynomial hierarchy. Later, Gaussian boson sampling (GBS), a variant of BS which avoids some of the costs of producing photon sources, was proposed by Hamilton, Kruse, Sansoni, Barkhofen, Silberhorn, and Jex [2]. GBS has gained popularity quickly, and various teams have made considerable efforts to build photonic quantum computers based on it [3]–[6].
The main proposed applications of GBS are to solve hard graph problems, such as counting perfect matchings [7], finding densest subgraphs [8], graph isomorphisms [9], planted bipartite cliques [10], and graph colourings [11]. The implementation in [4] uses samples from GBS for related graph search problems.
We start by defining the notation that we will use to sample from the GBS distribution on graphs (Theorem [thm:main-GBS]). While our presentation is self-contained, we mostly use the notation of [12]. We will subsequently extend to a related Boson sampling setting with non-negative weights (Theorem [thm:main-BS]). Given a graph \(G=(V,E)\) and a subset \(S \subseteq V\), let \(|S|\) denote the size of \(S\) and let \(PM(S)\) denote the number of perfect matchings in the induced subgraph \(G[S]\). There is a positive real number \(c\) — specifically \(c\) is the inverse of the maximum norm of an eigenvalue of the adjacency matrix of \(G\). Then the GBS distribution corresponding to \(G\) is defined as follows: \[\begin{align} \label{eqn:GBS-G-intro} \mu_{GBS,G}(S)\propto c^{2\left\vert S\right\vert}PM(S)^2. \end{align}\tag{1}\] In GBS applications the number \(c\) is in the range \((0,1)\) but in this paper we will not make any assumptions except that \(c\) is real and positive.
While there are certainly fast quantum algorithms for sampling from the GBS distribution — this is the whole point of GBS — it was an open question whether there is a fast classical algorithm, even in the unweighted (graph) case. Zhang, Zhou, Wang, Wang, Yang, Yang, Xue, and Li [12] showed recently that there is a fast algorithm for very dense graphs. In particular, they [12] showed that for any real number \(\xi\geq 1\) there is an algorithm which samples from \(\mu_{GBS,G}\) in time \(\widetilde{O}(n^{2\xi+18})\) for \(n\)-vertex graphs with minimum degree at least \(n-\xi\). Note that this is a very severe restriction on the density of the graph.
A related result which is interesting for graph applications is the quantum-inspired algorithm by Oh, Fefferman, Jiang, and Quesada [13]. Motivated by the application of GBS to graph problems and by the desire to develop efficient classical algorithms, they proposed sampling from a distribution that is a little different from \(\mu_{GBS,G}\), in the sense that the probability of the output \(S\) is proportional to \(PM(S)\) rather than to \(PM(S)^2\) as in 1 . They found that, in practice, their classical algorithm (from the different distribution) did not perform much worse in experiments than GBS. Their distribution can be sampled by two-photon boson sampling, which can be simulated classically (or can be directly sampled classically using the algorithm of Jerrum and Sinclair; see Proposition 2). Nevertheless, their work left open the question of whether there is an efficient classical algorithm for GBS in the non-negative case.
Our main result is that there is indeed a classical polynomial-time algorithm to sample from \(\mu_{GBS,G}\) (for all graphs).
theoremmainGBS There is an algorithm that, given a graph \(G=(V,E)\) with associated real \(c>0\), and arbitrary \(\varepsilon\), samples from a distribution that is \(\varepsilon\)-close to \(\mu_{GBS,G}\) in total variation distance, in time \(O(\overline{c}mn^4 \log^2 (n\overline{c}/\varepsilon))\) where \(m=\left\vert E\right\vert\), \(n=\left\vert V\right\vert\), and \(\overline{c}=\max\{1,c\}\).
Note that Theorem [thm:main-GBS] works for any \(c>0\), not just for \(c\in(0,1)\) as in the GBS application.
Remark 1. As a consequence of Theorem [thm:main-GBS], we have shown that there is no exponential-time quantum speedup for any application based on sampling from the distribution \(\mu_{GBS,G}\).
Our algorithm to sample from \(\mu_{GBS,G}\) is different from the approach in [12]. In fact, our algorithm is simpler, and faster, and it also applies to all graphs. We first construct the Cartesian product \(G\mathbin{\text{\scalebox{.84}{\square}}}K_2\) of \(G\) with an edge. We then note that the distribution \(\mu_{GBS,G}\) can be obtained by first sampling a perfect matching from \(G\mathbin{\text{\scalebox{.84}{\square}}}K_2\) with an appropriate weight function on edges and then projecting this perfect matching onto the vertex set matched by the original edges of \(G\). To sample weighted perfect matchings from \(G\mathbin{\text{\scalebox{.84}{\square}}}K_2\), we use the Markov chain of Jerrum and Sinclair [14]. The key property which guarantees that the Jerrum–Sinclair (JS) chain generates a perfect matching in polynomial time is that the ratio between the number of near-perfect matchings of \(G \mathbin{\text{\scalebox{.84}{\square}}}K_2\) (matchings with two unmatched vertices) and the number of perfect matchings of \(G\mathbin{\text{\scalebox{.84}{\square}}}K_2\) is bounded above by a polynomial in the size of the graph. We show that, for any graph \(G\), this property holds for \(G\mathbin{\text{\scalebox{.84}{\square}}}K_2\).
In addition to Gaussian boson sampling, we consider the original boson sampling model [1]. Here (using the notation of [15]), the input is an \(m\times n\) matrix \(A\) with \(n \leq m\), where \(A\) is the first \(n\) columns of a random \(m\times m\) unitary matrix. Let \(\Phi_{m,n}\) denote the set of non-negative integer vectors \(\mathbf{z}=\{z_1,...,z_m\}\) such that \(\sum_{i=1}^m z_i=n\). The boson sampling distribution is \[\begin{align} \label{eqn:BS-intro} \forall \mathbf{z}\in\Phi_{m,n},\quad\quad\mu_{BS}(\mathbf{z})=\frac{\left\vert\mathop{\mathrm{Perm}}(A_{\mathbf{z}})\right\vert^2}{\prod_{i=1}^m z_i!}, \end{align}\tag{2}\] where \(\mathop{\mathrm{Perm}}(\cdot)\) is the permanent function and \(A_{\mathbf{z}}\) is the square matrix composed of \(z_i\) copies of row \(i\) of \(A\). Exponential time classical simulation algorithms [15], [16] are known for this distribution, but it is conjectured that no polynomial-time classical approximate sampler exists [1]. Our second main result shows that, if the input is a non-negative matrix instead of a truncated unitary matrix, then one can efficiently sample from a similar distribution, where we modify the definition in 2 by renormalising appropriately, but retaining the fact that \(\forall \mathbf{z}\in\Phi_{m,n}\), \(\mu_{BS}(\mathbf{z})\propto\frac{\left\vert\mathop{\mathrm{Perm}}(A_{\mathbf{z}})\right\vert^2}{\prod_{i=1}^m z_i!}\).
theoremmainBS There is an algorithm that, given a non-negative real \(m\times n\) matrix \(A\) with \(n<m\) and a real number \(\varepsilon\in (0,1)\), samples from a distribution that is \(\varepsilon\)-close to \(\mu_{BS}\) in total variation distance, in time \(O\left(\frac{m^7n^{14}}{\varepsilon^7}\log^4\left(\frac{m n}{\varepsilon}\right)\right)\).
Our technique for proving Theorem [thm:main-BS] is a combinatorial reduction similar to the proof of Theorem [thm:main-GBS]. Given the matrix \(A\), we construct a bipartite graph \(G\) such that a certain distribution on weighted perfect matchings of \(G\) induces a distribution that is close to \(\mu_{BS}\). We then use the algorithm of Jerrum, Sinclair, and Vigoda [17] to efficiently sample from the distribution on perfect matchings.
When \(A\) is not a non-negative matrix, the Jerrum–Sinclair–Vigoda algorithm (JSV) is no longer applicable. However our construction would still work, if \(\left\vert\mathop{\mathrm{Perm}}(A_{\mathbf{z}})\right\vert\) could be efficiently approximated for each \(\mathbf{z}\). The task of approximating the norm of these permanents is thus the main source of difficulty for classically simulating boson sampling. The problem of approximating the permanent is NP-hard in general, even for real positive semi-definite matrices [18].
In Section 2 we define the terminology and review some known efficient classical algorithms for sampling matchings in graphs. We first prove Theorem [thm:main-BS] in Section 3 as a warmup, and then prove Theorem [thm:main-GBS] in Section 4.
2 Preliminaries↩︎
2.1 Classical algorithm to sample (perfect) matchings↩︎
Let \(G=(V,E)\) be a graph. A matching \(M\) of \(G\) is a subset of edges such that no two edges of \(M\) share an endpoint. Let \(V_M\) denote the vertex set of \(M\). If \(V_M=V\), then \(M\) is a perfect matching. If \(|V_M| = |V|-2\) then \(M\) is a near-perfect matching. Denote by \(\mathcal{M}\) the set of all matchings of \(G\), and by \(\mathcal{M}_k\) the set of matchings of size \(k\). Then, when \(\left\vert V\right\vert\) is even, the set of perfect matchings is \(\mathcal{M}_{\left\vert V\right\vert/2}\). Let \((\lambda_e)_{e\in E}\) be a collection of (non-negative) weights on the edges in \(E\). Define the following distribution over the matchings of \(G\): \[\begin{align} \label{eqn:matching} \forall M\in\mathcal{M},\quad \quad\mu_{matching,\lambda}(M)\propto \prod_{e\in M}\lambda_e. \end{align}\tag{3}\]
Jerrum and Sinclair [14] showed that there is a polynomial-time algorithm to approximately sample from the distribution 3 . To explain this more precisely, we need the following notation. The total variation (TV) distance between two distributions \(\mu\) and \(\nu\) over some discrete space \(\Omega\) is defined as \[\begin{align} \mathrm{dist}_{TV}(\mu,\nu):=\frac{1}{2}\sum_{\omega\in\Omega}\left\vert\mu(\omega)-\nu(\omega)\right\vert. \end{align}\]
Proposition 2. (Jerrum and Sinclair [14]) There is an algorithm that, given \(G=(V,E)\), weights \((\lambda_e)_{e\in E}\), and a real number \(\varepsilon\in (0,1)\), samples from a distribution that is \(\varepsilon\)-close to \(\mu_{matching,\lambda}\) in TV distance, in time \(O(\overline{\lambda}m n^2 \log \frac{n \overline{\lambda}}{\varepsilon})\), where \(m=\left\vert E\right\vert\), \(n=\left\vert V\right\vert\), and \(\overline{\lambda}=\max_{e\in E}\{1,\lambda_e\}\).
This particular running time is derived from [19].
Proposition 3 shows that the number of matchings of size \(k\) form a log-concave sequence (see [14]).
Proposition 3. For any graph \(G=(V,E)\) and \(1\le k\le \left\vert V\right\vert/2\), \(\left\vert M_{k-1}\right\vert\left\vert M_{k+1}\right\vert\le \left\vert M_k\right\vert^2\).
Later, to prove 12 , we will need a weighted version of Proposition 3, due to Heilmann and Lieb [20]. See Proposition 5.
Note that Proposition 2 samples weighted matchings but not perfect matchings. To efficiently sample perfect matchings, one can use the Jerrum-Sinclair algorithm only if \(\frac{\left\vert M_{t-1}\right\vert}{\left\vert M_{t}\right\vert}\) is bounded by a polynomial of the input size, where \(t=\left\vert V\right\vert/2\). When this is the case, due to log-concavity (Proposition 3), one can tune the edge weights so that perfect matchings show up sufficiently frequently. This will be our strategy for Theorem [thm:main-GBS] later.
On the other hand, in bipartite graphs, Jerrum, Sinclair, and Vigoda [17] showed that one can efficiently sample from the distribution 3 restricting to perfect matchings, without further assumptions. The running time of their algorithm was subsequently improved by Bezáková, Štefankovič, Vazirani, and Vigoda [21]. More precisely, for any graph \(G=(V,E)\) and any \(M \in \mathcal{M}_{\left\vert V\right\vert/2}\), let \[\begin{align} \label{eqn:PM-dist} \mu_{PM,\lambda}(M)\propto \prod_{e\in M}\lambda_e. \end{align}\tag{4}\]
Proposition 4. (Jerrum, Sinclair, Vigoda [17]; Bezáková, Štefankovič, Vazirani, and Vigoda [21]) There is an algorithm that, given a bipartite \(G=(V,E)\), weights \((\lambda_e)_{e\in E}\), and a real number \(\varepsilon\in (0,1)\), samples from a distribution that is \(\varepsilon\)-close to \(\mu_{PM,\lambda}\) in TV distance, in time \(O(n^7\log^4n+n^5\log\frac{n}{\varepsilon})\), where \(n=\left\vert V\right\vert\).
This running time comes from [21]. One needs to first spend \(O(n^7\log^4n)\) time to estimate certain weights, and then each subsequent sampling takes \(O(n^5\log\frac{n}{\varepsilon})\) time. See [21].
2.2 Permanent and Hafnian↩︎
Two matrix functions will be useful throughout this paper — the permanent and the hafnian. Let \(A\) be an \(n\times n\) matrix. The permanent is defined as follows \[\begin{align} \mathop{\mathrm{Perm}}(A):=\sum_{\sigma\in\mathcal{S}_{n}} \prod_{i=1}^n A_{i,\sigma(i)}, \end{align}\] where \(S_{n}\) is the symmetric group of order \(n\). If \(A\) is the biadjacency matrix of some bipartite graph \(G\), then \(\mathop{\mathrm{Perm}}(A)\) is the number of perfect matchings in \(G\).
Similarly, for a \(2n\times 2n\) symmetric matrix \(A\), \[\begin{align} \mathop{\mathrm{Haf}}(A) :=\frac{1}{2^n n!} \sum_{\sigma \in \mathcal{S}_{2n}} \prod_{i=1}^n A_{\sigma(2i-1),\sigma(2i)}. \end{align}\] If \(A\) is the adjacency matrix of some graph \(G\) with \(2n\) vertices, then \(\mathop{\mathrm{Haf}}(A)\) is the number of perfect matchings of \(G\).
3 Boson sampling↩︎
Boson sampling was proposed by Aaronson and Arkhipov [1] as a way to demonstrate quantum advantage. For exact classical simulation, see [15], [16]. Mathematically, the distribution to sample from is the following. Let \(A\) be an \(m\times n\) matrix with \(n\le m\). Let \(\mathbf{z}=\{z_1,\ldots,z_m\}\) be a vector of non-negative integers with \(\sum_{i=1}^{m} z_i=n\), and let \(\Phi_{m,n}\) denote the set of these vectors. Then, \[\begin{align} \label{eqn:BS} \forall \mathbf{z}\in \Phi_{m,n},\quad\quad\mu_{BS}(\mathbf{z})\propto \frac{\left\vert\mathop{\mathrm{Perm}}(A_\mathbf{z})\right\vert^2}{\prod_{i=1}^m z_i !}, \end{align}\tag{5}\] where \(A_{\mathbf{z}}\) is the square matrix composed of \(z_i\) copies of row \(i\) of \(A\). Physically, the matrix \(A\) is usually the first \(n\) columns of an \(m\times m\) Haar random unitary matrix, in which case the normalising factor in 5 is \(1\).
In the rest of this section, we prove Theorem [thm:main-BS], which we re-state for convenience.
Proof. Given the matrix \(A\), we will construct a bipartite graph \(G\) and weights \(\lambda\) such that its distribution \(\mu_{PM,\lambda}\) in 4 is \(\varepsilon/2\)-close to the distribution \(\mu_{BS}\) in 5 . We will then finish by invoking Proposition 4 with error \(\varepsilon/2\).
Let \(k = \lceil 4 n^2/\varepsilon\rceil\). We define \(G\) as follows: for each \(\ell \in \{1,2\}\) let \(L_\ell = \{ v_j^{(\ell)} : j\in [n]\}\). Let \(R_\ell = \{ u_{i,t}^{(\ell)}: i \in [m], t\in [k] \}\). The vertex set of \(G\) is \(L_1 \cup R_1 \cup L_2 \cup R_2\). The edge set of \(G\) is defined as follows: For each \(\ell\in \{1,2\}\) let \(E_\ell = \{(v_j^{(\ell)}, u_{i,t}^{(\ell)}): i\in [m], j\in [n], t\in [k], A_{i,j} \neq 0\}\). For each \(e = (v_j^{(\ell)}, u_{i,t}^{(\ell)})\in E_\ell\), \(\lambda_e = A_{i,j}\). Let \(E_R = \{(u_{i,t}^{(1)},u_{i,t}^{(2)}): j \in [m], t\in [k]\}\). For each \(e\in E_R\), \(\lambda_e = 1\). The edge set of \(G\) is \(E_1 \cup E_2 \cup E_R\). An example is given in Figure 1.
Let \(\mu_{PM,\lambda}\) from 4 be the distribution over perfect matchings in \(G\) with weights given by \(\lambda\). Let \(M\sim\mu_{PM,\lambda}\) be a sample from this distribution. Let \(S_1(M)\) be the set of vertices in \(R_1\) that are matched to vertices in \(L_1\) by \(M\) (rather than to vertices in \(R_2\)). Note that, because every vertex in \(L_1\) needs to be matched with some vertex in \(R_1\), \(\left\vert S_1(M)\right\vert=\left\vert L_1\right\vert=n\). We identify \(S_1(M)\) with a vector \(\mathbf{z}\in\Phi_{m,n}\), by defining \(z_i\) as follows for each \(i\in [m]\): \(z_i = \left\vert S_1(M) \cap\left\{u_{i,1}^{(1)},\ldots,u_{i,k}^{(1)}\right\}\right\vert\). Let \(\nu\) be the resulting distribution over \(\Phi_{m,n}\). We claim that, for \(k=\lceil \frac{4n^2}{\varepsilon}\rceil\), \[\begin{align} \label{eqn:closeness} \mathrm{dist}_{TV}(\nu,\mu_{BS}) \le \frac{\varepsilon}{2}. \end{align}\tag{6}\] The theorem follows from the claim and Proposition 4.
In the rest of the proof we show 6 . For any particular \(\mathbf{z}\in\Phi_{m,n}\), the way to obtain the sample \(\mathbf{z}\) is to first choose \(z_i\) copies of each \(u_i\) from \(\{u_{i,1}^{(1)},\ldots,u_{i,k}^{(1)}\}\) - denote the set of these by \(S_1\) - and then match \(S_1\) with \(L_1\) perfectly using edges in \(E_1\). Moreover, the vertices in \(R_1 \setminus S_1\) must be matched through \(E_R\) edges (which have weight \(1\)) to \(R_2\). As a result, the vertices in \(R_2\) with the same indices as those in \(S_1\) must be matched via edges in \(E_2\) to \(L_2\). Let \(S_2\) denote the set of these vertices in \(R_2\). The total weight of matching \(L_1\) with \(S_1\) perfectly is \(\mathop{\mathrm{Perm}}(A_\mathbf{z})\) since the rows of \(A_\mathbf{z}\) contain \(z_i\) copies of row \(i\) of \(A\). The matching from \(L_2\) to \(S_2\) is independent of this and also has total weight \(\mathop{\mathrm{Perm}}(A_\mathbf{z})\). To summarise, \[\begin{align} \nu(\mathbf{z})& \propto k^{-n} \nu(\mathbf{z})\propto k^{-n}\prod_{i=1}^m\binom{k}{z_i}\mathop{\mathrm{Perm}}(A_{\mathbf{z}})^2=\frac{\mathop{\mathrm{Perm}}(A_\mathbf{z})^2}{\prod_{i=1}^m z_i !}\cdot k^{-n}\prod_{i=1}^m\frac{k!}{(k-z_i)!}\notag\\ &= \frac{\mathop{\mathrm{Perm}}(A_\mathbf{z})^2}{\prod_{i=1}^m z_i !}\cdot k^{-n}\prod_{i=1}^m \prod_{j=0}^{z_i-1} (k-j) = \frac{\mathop{\mathrm{Perm}}(A_\mathbf{z})^2}{\prod_{i=1}^m z_i !}\cdot \prod_{i=1}^m\prod_{j=0}^{z_i-1}\left(1-\frac{j}{k}\right). \label{eqn:nu-dist} \end{align}\tag{7}\] Furthermore, using the facts that \(k=\lceil \frac{4n^2}{\varepsilon}\rceil\) and \(z_i \leq k\) and \(\sum_{i=1}^m z_i=n\), \[\begin{align} 1\ge \prod_{i=1}^m\prod_{j=0}^{z_i-1}\left(1-\frac{j}{k}\right) \ge \left(1-\frac{n}{k}\right)^n \ge e^{-2n^2/k} = e^{-\varepsilon/2}. \end{align}\] Plugging the estimate above into 7 , we have that the multiplicative error for each \(\mathbf{z}\) between the proportional weights of \(\nu\) and \(\mu_{BS}\) is between \(1\) and \(e^{-\varepsilon/2}\), which implies 6 . ◻
4 Gaussian boson sampling↩︎
Gaussian boson sampling (GBS) is a variant proposed by Hamilton, Kruse, Sansoni, Barkhofen, Silberhorn, and Jex [2]. See also [22]. For \(n\) photons and \(m\) modes, the input is a Gaussian state characterised by a \(2m\times 2m\) covariance matrix \(\sigma\). Let \(A\) be the sampling matrix given by \(A= \begin{pmatrix} 0 & I_m\\ I_m & 0 \end{pmatrix}\left(I_{2m}-\sigma_Q^{-1}\right)\), where \(I_m\) denotes the \(m\times m\) identity matrix and \(\sigma_Q:=\sigma + {I_{2m}}/{2}\). For any \(\mathbf{z}\in \Phi_{m,n}\), the probability of the output pattern \(\mathbf{z}\) is \[\begin{align} \label{eqn:GBS} \mu_{GBS}(\mathbf{z}):=\frac{\mathop{\mathrm{Haf}}(A_{\mathbf{z}})}{\sqrt{\det (\sigma_Q)}\prod_{i=1}^m z_i!}, \end{align}\tag{8}\] where \(A_{\mathbf{z}}\) is the square matrix corresponding to \(\mathbf{z}\). We will consider the case where \(n=O(\sqrt{m})\), where, with high probability, \(z_i\le 1\) (see [22]), and in this case the matrix \(A_{\mathbf{z}}\) is the principal submatrix of \(A\) formed by keeping rows and columns \(i\) and \(m+i\) for each \(i\in [m]\) with \(z_i=1\).1
Brádler, Dallaire-Demers, Rebentrost, Su, and Weedbrook [7] proposed using GBS to solve graph problems. They showed that one can choose the covariance matrix \(\sigma\) and \(c>0\) such that, for any graph \(G\) with adjacency matrix \(A\), the distribution in 8 can be interpreted as a distribution over \(S\subseteq V\): \[\begin{align} \label{eqn:GBS-G} \mu_{GBS,G}(S) \propto c^{2|S|} \mathop{\mathrm{Haf}}(A_S)^2, \end{align}\tag{9}\] where \(A_S\) is the principal submatrix of the adjacency matrix \(A\) corresponding to \(S\). Note that \(\mathop{\mathrm{Haf}}(A_S)\) counts the number of perfect matchings in the induced subgraph \(G[S]\). In fact, the distribution induced by GBS will have a parameter c in the range \((0,1)\).
In the rest of this section, we prove our main result, Theorem [thm:main-GBS], which shows that the distribution in 9 can be (classically) sampled in polynomial time for any graph \(G\) and \(c>0\).
The sampling algorithm that we use to prove Theorem [thm:main-GBS] takes the input graph \(G=(V,E)\) and constructs the graph \(G\mathbin{\text{\scalebox{.84}{\square}}}K_2\), which is the Cartesian product of \(G\) and an edge. In order to construct \(G\mathbin{\text{\scalebox{.84}{\square}}}K_2\), we start with a copy \(G'=(V',E')\) of \(G\). Then \(E_0 = \{ (v,v') : v\in V\}\) and \(G\mathbin{\text{\scalebox{.84}{\square}}}K_2 = (V\cup V', E \cup E' \cup E_0)\). See Figure 2 for an example. Moreover, for each edge \(e\in E\cup E'\), let \(\lambda_e=c^2\), and for each \(e\in E_0\), let \(\lambda_e=1\).
Suppose that \(M\) is a sample from the perfect matching distribution 4 on the graph \(G\mathbin{\text{\scalebox{.84}{\square}}}K_2\) with this weight function \(\lambda\). Let \(S_M\subset V\) be the set of vertices that are endpoints of some edge in \(M\cap E\). We have the following lemma, which shows that the induced distribution of \(S_M\) is given by Equation 9 .
Lemma 1. The induced distribution of \(S_M\) is exactly \(\mu_{GBS,G}(S)\).
Proof. For \(M\sim\mu_{PM,\lambda}\), as it is a perfect matching of \(G\mathbin{\text{\scalebox{.84}{\square}}}K_2\), each vertex \(v\) in \(V\) is either matched with an edge in \(M\cap E\) or \(M\cap E_0\). If \(v\) is matched in \(E_0\) then it is matched with its copy \(v'\). Thus the set of endpoints of edges in \(M\cap E'\), denoted \(S'\), is the set of copies of vertices in \(S_M\). See Figure 3 for an illustration.
To summarise, the probability of outputting a particular \(S\subset V\) is \[\begin{align} \ifthenelse{\isempty{M\sim\mu_{PM,\lambda}}} {\mathop{\mathrm{Pr}}\left[S_M=S\right]} {\mathop{\mathrm{Pr}}_{M\sim\mu_{PM,\lambda}}\left[S_M=S\right]} & \propto c^{2\left\vert S\right\vert} \mathop{\mathrm{Haf}}[A_S] \mathop{\mathrm{Haf}}[A_{S'}]\\ & = c^{2\left\vert S\right\vert} \mathop{\mathrm{Haf}}[A_S]^2 \qedhere \end{align}\] ◻
Our problem now is to sample from \(\mu_{PM,\lambda}\) for \(G \mathbin{\text{\scalebox{.84}{\square}}}K_2\). For this, we use Proposition 2, which can be used to efficiently sample perfect matchings in only if the ratio between the number of perfect and near-perfect matchings is bounded from above by a polynomial. We will show that this is indeed the case for \(G\mathbin{\text{\scalebox{.84}{\square}}}K_2\).
To be more specific, define a new weight function \(\lambda'\) for \(G \mathbin{\text{\scalebox{.84}{\square}}}K_2\) as \(\lambda_e'=4n^2\lambda_e\). This weight function favours matchings of larger sizes, and will eventually be the weight function that we will use when invoking the Jerrum–Sinclair chain [14]. For a matching \(M\) of \(G\mathbin{\text{\scalebox{.84}{\square}}}K_2\), let \(w(M):=\prod_{e\in M}\lambda_e\) and \(w'(M):=\prod_{e\in M}\lambda'_e\). Let \(M_k\) denote the set of matchings with \(k\) edges in \(G\mathbin{\text{\scalebox{.84}{\square}}}K_2\), where \(0\le k\le n\). Let \(Z_k\) be the total weight from \(k\)-edge matchings under \(\lambda\), namely, \[\begin{align} Z_k:=\sum_{M\in M_k}w(M), \end{align}\] and \(Z:=\sum_{k=1}^n Z_k\) be the partition function of \(\mu_{matching,\lambda}\). Similarly, define \(Z_k'\) and \(Z'\) for \(\lambda'\). Note that the proof of Lemma 1 implies that \[\begin{align} \label{eqn:sum-Z95n} Z_{n} = \sum_{S\subseteq V}c^{2\left\vert S\right\vert}\mathop{\mathrm{Haf}}(A_S)^2. \end{align}\tag{10}\]
Lemma 2. \(Z_{n-1}< 2n^2 Z_{n}\).
Proof. We write \(\mathop{\mathrm{Haf}}(S)\) for \(\mathop{\mathrm{Haf}}(A_S)\). Note that by choosing the two vertices to remove, we have \[\begin{align} Z_{n-1} &= 2 \sum_{(u_1,u_2) \in \binom{V}{2}} \sum_{S \subseteq V \setminus \{u_1,u_2\}} c^{2\left\vert S\right\vert+2}\mathop{\mathrm{Haf}}(S) \cdot \mathop{\mathrm{Haf}}(S \cup \{u_1,u_2\})\notag\\ &\quad +\sum_{u_1,u_2 \in V} \sum_{S \subseteq V \setminus \{u_1,u_2\}} c^{2\left\vert S\right\vert+2}\mathop{\mathrm{Haf}}(S \cup \{u_1\}) \cdot \mathop{\mathrm{Haf}}(S \cup \{u_2\}).\label{eqn:Z95n-1-sum} \end{align}\tag{11}\] By the AM-GM inequality, it holds that \[\begin{align} &\quad \sum_{(u_1,u_2) \in \binom{V}{2}} \sum_{S \subseteq V \setminus \{u_1,u_2\}} c^{2\left\vert S\right\vert+2}\mathop{\mathrm{Haf}}(S) \cdot \mathop{\mathrm{Haf}}({S \cup \{u_1,u_2\}}) \\ &\le \frac{1}{2}\sum_{(u_1,u_2) \in \binom{V}{2}} \sum_{S \subseteq V \setminus \{u_1,u_2\}} \left(c^{2\left\vert S\right\vert}\mathop{\mathrm{Haf}}(S)^2+ c^{2\left\vert S\right\vert+4}\mathop{\mathrm{Haf}}(S \cup \{u_1,u_2\})^2\right)\\ &\le \binom{n}{2} \sum_{S \subseteq V} c^{2\left\vert S\right\vert}\mathop{\mathrm{Haf}}(S)^2 = \binom{n}{2}Z_{n}, \end{align}\] where we used 10 in the last line. Similarly, \[\begin{align} &\quad \sum_{u_1,u_2 \in V} \sum_{S \subseteq V \setminus \{u_1,u_2\}} c^{2\left\vert S\right\vert+2}\mathop{\mathrm{Haf}}(S \cup \{u_1\}) \cdot \mathop{\mathrm{Haf}}(S \cup \{u_2\}) \\ & \le \frac{1}{2} \sum_{u_1,u_2 \in V} \sum_{S \subseteq V \setminus \{u_1,u_2\}} \left(c^{2\left\vert S\right\vert+2}\mathop{\mathrm{Haf}}(S \cup \{u_1\})^2+ c^{2\left\vert S\right\vert+2}\mathop{\mathrm{Haf}}(S \cup \{u_2\})^2\right)\\ &\le n^2 \cdot \sum_{S \subseteq V} c^{2\left\vert S\right\vert}\mathop{\mathrm{Haf}}(S)^2=n^2 Z_{n}. \end{align}\] Combining with 11 , it follows that \[\begin{align} Z_{n-1} &< 2 n^2 Z_{n}. \qedhere \end{align}\] ◻
In addition, Heilmann and Lieb [20] showed that the sequence \((Z_k)\) is log-concave. This is a weighted version of Proposition 3.
Proposition 5. For any \(1\le k\le n-1\), \(Z_{k-1}Z_{k+1}\le Z_{k}^2\).
By Proposition 5 and Lemma 2, we have that \(Z_{k-1}\le 2n^2Z_k\) for any \(1\le k\le n\). Since \(Z_k'=(4n^2)^k Z_k\), it follows that for any \(1\le k\le n\), \(Z_{k-1}'\le Z_k'/2\). Thus, \[\begin{align} \label{eqn:PM-percent} Z'=\sum_{i=0}^nZ_i'\le \sum_{i=0}^n 2^{i-n}Z_n' <2Z_n'. \end{align}\tag{12}\]
We can now prove Theorem [thm:main-GBS].
Proof. By Lemma 1, we just need to sample from \(\mu_{PM,\lambda}\) over perfect matchings in \(G\mathbin{\text{\scalebox{.84}{\square}}}K_2\). To this end, we use Proposition 2 to approximately sample from the distribution \(\mu_{matching,\lambda'}\) in 3 for \(G\mathbin{\text{\scalebox{.84}{\square}}}K_2\) with the weight function \(\lambda'\) as described above. We set the error as \(\varepsilon'=\min\{1/4,\varepsilon/2\}\) in Proposition 2. By 12 , there is \(1/2-\varepsilon'\geq1/4\) probability that the algorithm outputs a perfect matching. We keep rejecting until we see a perfect matching. The distribution \(\mu_{matching,\lambda'}\) conditioned on outputting a perfect matching is exactly \(\mu_{PM,\lambda}\). Thus we can approximately sample from \(\mu_{PM,\lambda}\) with \(\varepsilon\) error after at most \(O(\log{(1/\varepsilon)})\) rejections. As for the running time, we plug in \(\lambda = O(n^2)\) in Proposition 2, which finishes the proof. ◻
5 Concluding remarks↩︎
In this paper we showed that there are polynomial-time classical algorithms to (approximately) sample from the distribution 9 , related to Gaussian Boson Sampling on graphs, and from the distribution 5 when the input matrix \(A\) is non-negative. It would be interesting to know whether the non-negative restriction can be relaxed. While the algorithm on which our method is based works only for non-negative matrices, in the more general case, our construction actually enables a reduction from sampling to approximating permanents of some complex weighted matrices. This is because, as can be seen from the proof of Theorem [thm:main-BS], the marginal probability of choosing a row is the ratio between the permanents of two complex weighted matrices, even when conditioned on previous choices. Thus, one can sample from the (conditional) marginal distributions of the rows, one by one, using an oracle that approximates the permanent of complex weighted matrices. However, the latter task is NP-hard in general [18], and it is an interesting open question to what extent this strategy can be useful.
Acknowledgement↩︎
HG would like to thank Raul Garcia-Patron for pointing out a few useful references.
This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 947778).