{"id":"4eec46ea-25e7-44ce-a8e3-b577f74477bc","arxiv_id":"1908.08068","paper_version":3,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":7.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":2,"one_line_summary":"A chain-rule algorithm samples each mode of a Gaussian Boson Sampler sequentially, giving exact simulation in polynomial space and time exponential in the detected photon number.","lead":"This paper introduces an exact classical algorithm for simulating Gaussian Boson Sampling with polynomial memory and time that grows exponentially with the detected photon number. It provides a practical benchmark for testing photonic quantum devices and extends to threshold detectors, displaced Gaussian states, and mixtures.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"The algorithm as written is not exact: it truncates photon counts at nmax, and the complexity analysis omits the growth of nmax required by the error bound in Appendix B.","rationale":"The paper's main contribution is an 'exact' classical sampler with polynomial space and time exponential in photon number. For that claim to hold, the sampler must sample from the exact GBS distribution. The algorithm as written in Section II.B does not: it truncates at nmax and samples from the normalized truncated distribution. The Appendix B bound only gives a total-variation guarantee for the truncated distribution; it does not make the algorithm exact. Moreover, the complexity analysis treats nmax as a constant, but for a fixed error or for exactness, nmax has to grow with m (and possibly with N) because the bound in Eq. (B2) sums m mode tails. This means the O(m N^3 2^N) time bound is only valid if one accepts a distribution that differs from the ideal GBS distribution by a controlled but nonzero error, and the constant depends on the error and on m. The reader's verdict already flags this as the weakest assumption, and my analysis agrees with that identification. The secondary concern about the Appendix A non-negative kernel proof (the apparent circularity in identifying O=N) is real but less central because it affects only the approximate polynomial-time variant, not the main exponential-time sampler. I do not see a fatal flaw in the chain-rule/hafnian construction itself; with a proper inverse-CDF strategy and an analysis of expected runtime, the approach could likely be made exact in a well-defined sense. The appropriate verdict therefore remains CONDITIONAL: the authors should revise the claims to distinguish exact sampling (with a rigorous unbounded-support procedure) from the implemented truncated approximation, and the complexity analysis should account for the dependence of nmax on problem size and target error. The changes would not overturn the value of the contribution, but they are necessary before the 'exact' label is accepted.","tokens_in":12004,"tokens_out":6962,"duration_ms":70210,"concrete_test":"Analytical check: using Eq. (B2), compute the minimum nmax needed for total-variation error epsilon=0.01 for m=100 and m=1000 modes at unit mean photon number per mode. If nmax grows with m, then 2^{nmax} is m-dependent and the O(m N^3 2^N) form in Eq. (13) is incomplete. Complement this by running the supplied Walrus sampler on a single squeezed mode and comparing the empirical photon-number histogram to the exact thermal distribution for several nmax values; a nonzero total-variation distance directly shows the implementation samples a truncated distribution.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The central claim of exact simulation is not supported by the algorithm in Section II.B. Step 1 (and analogously steps 2-3) computes p(s_k=0), ..., p(s_k=nmax) and samples from the resulting truncated distribution. This is an approximation, not exact sampling from the unbounded GBS distribution. Appendix B acknowledges this by bounding the total-variation distance of the truncated distribution; it does not establish exactness. The complexity analysis in Section II.C (Eqs. 12-13) treats nmax as a fixed constant, absorbing 2^{nmax} into the O-constant. However, for a fixed error epsilon, the bound in Eq. (B2) forces nmax to grow with the number of modes m; for squeezed-state marginals the tail decays exponentially, so nmax = O(log(m/epsilon)). Then the hafnian cost in Eq. (12) carries a factor 2^{nmax} = O(m/epsilon), and the claimed O(m N^3 2^N) time does not hold for an exact sampler; the dependence on m is hidden in the constant. The paper does not present an alternative inverse-CDF procedure that samples from the full unbounded support without a cutoff, nor does it analyze such a procedure's runtime. The 'exact' label therefore applies only to a truncated distribution, and the stated complexity omits the cost of controlling truncation error.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper presents a classical sampling algorithm for Gaussian Boson Sampling (GBS). The method samples modes sequentially using the chain rule of probability, computing each conditional probability from the reduced covariance matrix of the modes already sampled and the corresponding hafnian. The authors claim exact sampling in O(m^2) space and O(m N^3 2^N) time, where N is the number of detected photons, and argue that for fixed mode number the cost is equivalent to that of computing output probabilities up to constant factors. They extend the method to threshold detectors, displaced GBS, linear combinations of Gaussian states, and a polynomial-time approximate variant for non-negative kernel matrices. The paper also reports an implementation in The Walrus library and benchmarks it on Haar-random interferometers and Erdős–Rényi graph encodings.","tokens_in":12289,"tokens_out":5731,"duration_ms":63160,"significance":"If the central claims were fully established, the paper would be a useful contribution: it would show that classical simulation of GBS has space complexity only quadratic in the mode number and time complexity governed by the cost of hafnian computation, placing GBS simulation on a footing comparable to Boson Sampling. The paper also provides a publicly available implementation, benchmarks with concrete runtimes, and a clear numerical demonstration of the challenges facing the approximate non-negative-kernel variant. These strengths make the work potentially valuable to the quantum-simulation and quantum-advantage community. However, the exactness and complexity claims are not supported as stated because the implemented algorithm truncates the photon-number distribution, and the Appendix A proof underlying the polynomial-time variant contains a logical gap.","major_comments":[{"comment":"The algorithm as described is not an exact sampler from the GBS distribution. In Step 1, and analogously in Steps 2 and 3, the algorithm computes p(s_k = 0), ..., p(s_k = nmax) and samples from the resulting truncated distribution, rather than from the full unbounded photon-number distribution. Equation (11) is exact only if every conditional distribution is sampled over its full support. Appendix B explicitly acknowledges this by bounding the total-variation distance of the truncated distribution, but a bound on the approximation error does not establish exactness. The manuscript should either provide a procedure that samples exactly from the unbounded support or revise the title, abstract, and Section II.B to describe an approximate sampler with a controlled error, and should state precisely how the remaining probability mass (the event mapped to ∅ in Appendix B) is treated in the sampling step.","section":"II.B, Steps 1–3; Appendix B"},{"comment":"The complexity analysis treats nmax as a fixed constant and absorbs 2^{nmax} into the O-constant, but this is not legitimate for the stated claims. If the sampler is required to approximate the GBS distribution within a fixed total-variation error, Eq. (B2) forces nmax to grow with the number of modes m (and with the state parameters), so the per-mode cost contains a factor 2^{nmax} that depends on m. Consequently, the claimed O(m N^3 2^N) complexity for an exact or fixed-error sampler is not established. The paper should state the complexity as a function of nmax and the target error, or restrict the exactness/complexity claim to a regime where m is fixed and nmax is chosen accordingly, with that dependence made explicit.","section":"II.C, Eqs. (12)–(13)"},{"comment":"The proof that a non-negative proper covariance matrix Q leads to a non-negative kernel matrix A contains a logical gap. The identification Q^{-1} = 1 - O and the comparison with the M-matrix decomposition Q^{-1} = λ0 1 - N, setting λ0 = 1, assumes part of what is to be proved. The chain of inequalities in Eq. (A9) is in the Loewner (positive-semidefinite) order, and positive-definiteness bounds do not imply entrywise non-negativity of O. Since the polynomial-time approximate algorithm of Section II.D.4 depends on this lemma, the proof needs to be replaced by a valid argument or a precise reference establishing the claim.","section":"Appendix A, Eqs. (A3)–(A9)"}],"minor_comments":[{"comment":"The reported average error ε = 63% and maximum error ε = 228% in the approximate non-negative-kernel algorithm are very large; the text acknowledges this, but the figure caption and the phrase 'approximate sampling algorithm' should more prominently warn that the sampled distribution can be far from the target GBS distribution at the tested parameters.","section":"III.B, Fig. 2 and text"},{"comment":"The statement that 'for all values of N for which hafnians can be computed in practice, it is preferable to employ the algorithm of Ref. [38]' would benefit from a quantitative comparison or a reference to the actual crossover point, since it is not obvious from the asymptotic exponents alone.","section":"II.C, paragraph after Eq. (13)"},{"comment":"The notation 'S > nmax' is introduced informally; a short formal definition (e.g., S in N^m with at least one coordinate exceeding nmax) would improve readability.","section":"Appendix B, Eq. (B2)"}],"recommendation":"major_revision","confidential_remarks":"The manuscript's key sampling derivation via the chain rule and reduced covariance matrices is sound, and the implementation is a useful practical contribution. The main obstacle is the mismatch between the claimed 'exact' simulation and the truncated sampler actually described; this is fixable by reframing the paper around an approximate sampler with explicit error bounds and a revised complexity expression. I would also ask an expert in matrix analysis to verify the Appendix A lemma, since the current proof appears to rely on an unjustified identification. Some key formulas are cited to the authors' own previous work, which is not unusual, but independent verification of the non-negativity lemma would strengthen the manuscript."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Colleague,\n\nThe thing to know: this is the first GBS sampling algorithm with photon-number-resolving detectors that runs in polynomial memory and time exponential in photon number, and the core chain-rule idea is correct. Worth engaging.\n\nWhat's new: Clifford and Clifford gave the chain-rule approach for standard boson sampling; the authors adapt it to GBS using Gaussian marginals and hafnians. The reduced covariance matrix trick gives O(m^2) space and O(m N^3 2^N) time for fixed mode count and fixed photon cutoff. That is a real improvement over the threshold-only sampler, which needed exponential memory. They ship code in The Walrus and benchmark both the exact-ish and approximate variants.\n\nWhere it gets soft: the word 'exact' is doing too much work. Section II.B explicitly truncates at nmax and samples from the resulting distribution; Appendix B only bounds the total variation distance to the true distribution. That makes the implemented algorithm an approximation, not exact sampling. The complexity analysis treats nmax as a constant, but the error bound forces nmax to grow roughly like log(m/epsilon) for squeezed states, so the stated O(m N^3 2^N) hides a factor that grows with m and 1/epsilon. The paper should either present an unbounded inverse-CDF sampler or reclassify the algorithm as approximate with a rigorous runtime that accounts for nmax. This is a fixable issue, not a fatal one: the chain rule itself is exact, and the truncation error can be controlled.\n\nThe reader flagged a circular step in Appendix A's non-negative kernel proof. I read the proof carefully and think the concern does not land: the argument that O=1-Q^{-1} has eigenvalues in (-1,1) follows from Q^{-1} positive definite with eigenvalues in (0,2), so lambda0=1 is justified. The decomposition Q^{-1}=lambda0 1 - N is not circular; it is the standard M-matrix representation.\n\nOne smaller caveat: the approximate polynomial-time algorithm for non-negative kernels has large errors in the benchmark (63% average) and is of limited practical use, but the paper reports this honestly.\n\nBottom line: for GBS practitioners and complexity theorists, this is a useful result and the first direct simulation method in polynomial space. Fix the exactness language and the nmax dependence, then it is solid. I'd send it to a serious referee.","headline":"First polynomial-memory GBS sampler with PNR detectors; the chain-rule construction is sound, but the 'exact' label is an overstatement once the finite photon cutoff is in play.","tokens_in":12799,"tokens_out":2659,"would_cite":true,"duration_ms":25630,"reading_group":"yes","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":[],"pacs":[],"model":"deepseek-v4-flash","headline":"Classical algorithm samples GBS exactly in polynomial space","keywords":["Gaussian Boson Sampling","classical simulation","hafnian","chain rule","conditional sampling","covariance matrix","threshold detectors","non-negative kernel matrix"],"falsifier":"Run the algorithm on a small GBS instance, say four modes with fixed squeezing and a known interferometer, where exact output probabilities can be enumerated by brute force, and increase the maximum photon number $n_{\\max}$. The central claim requires the total variation distance between the sampled and exact distributions to go to zero as $n_{\\max}$ grows; if the distance saturates above zero, or if the $n_{\\max}$ needed for a fixed accuracy forces a runtime larger than $O(mN^3 2^N)$ as the mode count grows, the exactness claim fails.","tokens_in":11782,"feed_emoji":"🎲","tokens_out":11486,"duration_ms":102251,"temperature":0.7,"pith_summary":"The paper introduces a classical algorithm that samples exactly from the distribution produced by Gaussian Boson Sampling (GBS), in which a multi-mode Gaussian state is measured with photon-number detectors. The method applies the chain rule of probability, sampling each mode's photon number conditional on the previously sampled numbers, and it computes each conditional probability from the reduced covariance matrix of the modes already processed. For $m$ modes and $N$ detected photons the algorithm runs in $O(m^2)$ space and $O(mN^3 2^N)$ time, so for a fixed number of modes it is, up to constants, as costly as computing a single GBS output probability. The same conditioning strategy covers GBS with threshold detectors, displacements, and linear combinations of Gaussian states, and it yields a polynomial-time approximate sampler when the kernel matrix is non-negative. This matters because claims that GBS is hard to simulate classically are asymptotic; the algorithm gives a concrete benchmark of what classical simulation actually costs.","feed_headline":"Classical algorithm samples GBS exactly in polynomial space","feed_subtitle":"For a fixed mode count, drawing one sample costs about the same as computing one output probability.","key_machinery":"The load-bearing object is the reduced covariance matrix $Q^{(k)}=V^{(k)}+1/2$ and the kernel matrix $A^{(k)}=X(1-(Q^{(k)})^{-1})$ for the first $k$ modes. The chain rule of probability, $p(s_1,\\dots,s_m)=\\prod_{k=1}^m p(s_k|s_1,\\dots,s_{k-1})$, is the engine: each conditional probability is a ratio of two marginal probabilities of the form in Eq. (8), and each marginal is a hafnian whose matrix is obtained from $A^{(k)}$ by repeating rows and columns according to the photon pattern. The hafnian, the sum over perfect matchings that counts the terms in the output probability, is the only exponential ingredient, since Gaussian conditioning keeps all reduced matrices polynomial in size. For the non-negative case, the hafnian is replaced by the Monte Carlo identity $\\mathrm{Haf}(A)=\\mathbb{E}[\\det W]$ for a random skew-symmetric $W$, which supplies the polynomial-time approximate variant.","core_discovery":"The central claim is that sampling from the GBS distribution reduces to a sequence of conditional probability computations. Writing $S^{(k)}=(s_1,\\dots,s_k)$ for a partial pattern, the paper uses $p(S^{(k)})=\\mathrm{Haf}(A_s^{(k)})/(\\sqrt{\\det Q^{(k)}}\\,s_1!\\cdots s_k!)$, where $Q^{(k)}=V^{(k)}+1/2$ is built from the reduced covariance matrix $V^{(k)}$ of the first $k$ modes and $A^{(k)}=X(1-(Q^{(k)})^{-1})$ is the corresponding kernel matrix; the hafnian $\\mathrm{Haf}$ is the sum over perfect matchings that gives the output probability. Because Gaussian marginals are Gaussian, every reduced matrix is obtained efficiently, and the algorithm samples each $s_k$ from $p(s_k|s_1,\\dots,s_{k-1})$ by computing two such marginal quantities. The paper argues this gives exact sampling with $O(m^2)$ space and $O(mN^3 2^N)$ time, where $N$ is the total number of detected photons, and extends the same chain-rule construction to threshold detectors, displaced Gaussian states, and linear combinations of Gaussian states. For non-negative kernel matrices, the hafnian is replaced by a randomized estimator of the form $\\mathbb{E}[\\det W]$, giving an approximate polynomial-time sampler.","pith_inferences":["Because the only property the sampler needs is efficient marginals, the method should extend to any state family with efficiently computable reduced descriptions, not just Gaussian states; the paper notes the Gaussian property but does not pursue this general criterion.","The runtime is driven by the number of detected photons, so the hardest instances are high-brightness outputs; low-photon samples, even in many modes, should be classically cheap to produce.","The non-negative case's practicality is gated by the $1/\\sqrt{M}$ Monte Carlo convergence of the hafnian estimator; a variance-reduced estimator or importance-sampling scheme would be a natural development, and the paper's own error benchmarks show why it is needed.","One can compute the Appendix B bound in polynomial time, so an instance-specific choice of maximum photon number could certify a target total-variation error; the tradeoff between that cutoff, the mode count, and the runtime is a direct next calculation."],"forward_implications":["For a fixed number of modes, exact classical simulation of GBS costs the same order as computing a single output probability, up to constant factors.","The $O(m^2)$ space requirement removes the exponential memory bottleneck of earlier threshold-detector simulation algorithms.","The same chain-rule sampler handles GBS with threshold detectors, displaced Gaussian states, and linear combinations of Gaussian states, with the time cost changed only by known prefactors.","When the kernel matrix is non-negative, approximate sampling runs in polynomial time; the paper's benchmark reports errors up to 228% in the normalization proxy, so the polynomial speedup carries a substantial accuracy cost."],"supporting_citations":[{"why":"Supplies the GBS output probability formula linking the kernel matrix to the hafnian, the foundation of the sampling algorithm.","marker":"[10]"},{"why":"Provides the $O(N^3 2^N)$ hafnian and loop-hafnian algorithms that set the per-photon time complexity.","marker":"[37]"},{"why":"Gives the hafnian algorithm with better constants that the implementation actually uses for practical sizes.","marker":"[38]"},{"why":"Introduces the torontonian for threshold-detector GBS, which the extension to threshold detectors relies on.","marker":"[26]"},{"why":"Supplies the polynomial-time randomized estimator for hafnians of non-negative matrices that powers the approximate algorithm.","marker":"[41, 42]"},{"why":"States the physical condition on the covariance matrix used in the Appendix A argument that non-negative covariances give non-negative kernels.","marker":"[45]"},{"why":"Provides the M-matrix property used to conclude the inverse covariance matrix has the sign pattern needed for the non-negative kernel claim.","marker":"[46]"}],"fun_headline_variants":["Exact GBS sampling in polynomial space, exponential time","GBS exact simulation: polynomial space, exponential time","Exact GBS samples with just polynomial memory","Sampling GBS exactly in polynomial space","One probability per sample: exact GBS in poly space"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The algorithm is exact only if the photon number in every mode is drawn from its full, unbounded range of possibilities; any practical maximum photon number turns the simulation into an approximation, and the stated runtime treats that maximum as a fixed constant rather than as a quantity that must grow with the number of modes.","fun_headline_variants_meta":{"raw":{"variants":["Exact GBS sampling in polynomial space, exponential time","GBS exact simulation: polynomial space, exponential time","Exact GBS samples with just polynomial memory","Sampling GBS exactly in polynomial space","One probability per sample: exact GBS in poly space"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000466,"raw_usage":{"total_tokens":2343,"prompt_tokens":982,"completion_tokens":1361,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":598,"completion_tokens_details":{"reasoning_tokens":1286}},"tokens_in":598,"tokens_out":1361,"duration_ms":9429,"temperature":1.0,"reasoning_tokens":1286,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-14T11:52:48.223982+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Run the algorithm on a small GBS instance, say four modes with fixed squeezing and a known interferometer, where exact output probabilities can be enumerated by brute force, and increase the maximum photon number $n_{\\max}$. The central claim requires the total variation distance between the sampled and exact distributions to go to zero as $n_{\\max}$ grows; if the distance saturates above zero, or if the $n_{\\max}$ needed for a fixed accuracy forces a runtime larger than $O(mN^3 2^N)$ as the mode count grows, the exactness claim fails.","supporting_citations":[{"cited_title":"Tillmann, B","cited_arxiv_id":null,"evidence_quote":"Supplies the GBS output probability formula linking the kernel matrix to the hafnian, the foundation of the sampling algorithm."},{"cited_title":"Weedbrook, S","cited_arxiv_id":null,"evidence_quote":"Provides the $O(N^3 2^N)$ hafnian and loop-hafnian algorithms that set the per-photon time complexity."},{"cited_title":"Point Processes with Gaussian Boson Sampling","cited_arxiv_id":"1906.11972","evidence_quote":"Introduces the torontonian for threshold-detector GBS, which the extension to threshold detectors relies on."},{"cited_title":"Bj¨ orklund, B","cited_arxiv_id":null,"evidence_quote":"States the physical condition on the covariance matrix used in the Appendix A argument that non-negative covariances give non-negative kernels."},{"cited_title":"Bj¨ orklund, inProceedings of the twenty-third annual ACM-SIAM symposium on Discrete Algorithms (Soci- ety for Industrial and Applied Mathematics, 2012) pp","cited_arxiv_id":null,"evidence_quote":"Provides the M-matrix property used to conclude the inverse covariance matrix has the sign pattern needed for the non-negative kernel claim."}],"review_version":1}