Derivation of the Loop Hafnian Generating Function for Arbitrary Symmetric Matrices via Gaussian Integration

Sergey V. Tarasov

Derivation of the Loop Hafnian Generating Function for Arbitrary Symmetric Matrices via Gaussian Integration

T0 review · 0 major / 5 minor · reviewed 2026-08-06 · deepseek-v4-flash

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗

Pith's one-line read For every symmetric complex matrix, the loop hafnian generating function previously tied to quantum Gaussian states follows from Gaussian integration alone.

desk verdict A short, clean proof that removes a known restriction on the loop hafnian generating function; worth a serious referee but not a paradigm shift. read the letter →

arxiv 2507.16100 v1 pith:QZCSI6MG submitted 2025-07-21 quant-ph cond-mat.stat-mech

Sergey V. Tarasov This is my paper

classification quant-phcond-mat.stat-mech MSC 15A1505A15

keywords loophafniangeneratingfunctionsymmetriccomplexmatrixGaussianintegrationmastertheoremquantumstatesmoments

verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

Loop hafnians are polynomial sums over pairings and singletons of a symmetric matrix's entries; they appear in counting and quantum-optics problems and are hard to compute in general. This paper proves that the exponential expression in Eq. (4) is their generating function for any even-size symmetric complex matrix and any complex vector, not just for matrices descending from quantum Gaussian states. The earlier restricted version came from photon-number probability sums; the proof here is Gaussian integration plus a derivative identity. If the theorem is right, the formula can be used broadly in simulations and combinatorial calculations without checking quantum-state symmetries.

What carries the argument

The load-bearing machinery is the shifted Gaussian-integral representation $$g(x)=\exp\left(\tfrac{1}{2}$x^{{T}}$Sx+$v^{{T}}$x\right)=\exp\left(-\tfrac{1}{2}$v^{{T}}$$S^{{-1}}$v\right)(2\pi)^{-m}\int_{\mathbb{R}^{2m}}\exp\left(-\tfrac{1}{2}\$xi^{{T}}$\xi+$v^{{T}}$$S^{{-1/2}}$\xi+$x^{{T}}$$S^{{1/2}}$\xi\right)$d^{{2m}}$\xi,$$ combined with the derivative lemma (Eq. (5)) that expresses every loop hafnian as mixed derivatives of $g$ at $x=0$. Inserting the block counter-diagonal matrix $Z=\begin{pmatrix}0&\operatorname{diag}\{z_j\}\\ \operatorname{diag}\{z_j\}&0\end{pmatrix}$ converts those derivatives into derivatives with respect to generating variables $z_j$, so the whole expression becomes an evaluable Gaussian integral. Continuity of both sides near $z=0$ then carries the identity from invertible diagonalizable $S$ to all symmetric $2m\times 2m$ matrices.

What would settle it

Take the $2\times 2$ non-diagonalizable complex symmetric matrix $S=\begin{pmatrix} i & 1\\ 1 & -i\end{pmatrix}$, a generic vector $v$, and expand the left side of Eq. (4) through degree two in $z_1$; compute the corresponding loop hafnians of the extended matrix directly from Definition 2.2. A mismatch between the two coefficient lists would refute the theorem, while agreement at this singular example would support the extension beyond the density argument.

Watch

Extended reading notes

Core claim

The central claim is the loop hafnian master theorem, Eq. (4): for every symmetric $2m\times 2m$ complex matrix $S=S^{T}$ and every $2m$-dimensional complex vector $v$, with $Z=\begin{pmatrix}0&\operatorname{diag}\{z_j\}\\ \operatorname{diag}\{z_j\}&0\end{pmatrix}$, the identity $$\frac{\exp\left(\frac{1}{2}$v^{{T}}$(1-ZS)^{-1}Zv\right)}{\sqrt{\det(1-ZS)}}=\sum_{\{n_j\}}\operatorname{lhaf}\left(\tilde{S}_{\{n_j\},\{n_j\}};\tilde{v}_{\{n_j\},\{n_j\}}\right)\prod_{j=1}^{m}\frac{$z_j^{{n_j}}$}{n_j!}$$ holds as an exponential generating function in the variables $z_j$. The earlier quantum-optical derivation was limited to matrices with the block symmetries of Gaussian-state covariance matrices; here that restriction is shown to be unnecessary. The proof replaces the quantum specificity by a shifted Gaussian integral and ends with a continuity argument that reaches every symmetric even-size complex matrix.

Load-bearing premise

The proof writes down the Gaussian integral under the assumption that $S$ is invertible and diagonalizable with a symmetric square root; the step from there to every symmetric matrix relies on continuity of both sides near $z=0$.

Editorial extensions

If this is right

Any use of the generating function in algorithms that simulate lossy or partially distinguishable optical circuits is now justified for arbitrary symmetric matrices, removing the need to verify Gaussian-state block symmetries or to treat the formula as a numerically supported assumption.
Loop hafnians of odd-size matrices are covered by the same generating function through the relation $\operatorname{lhaf} S_{2n-1}=\operatorname{lhaf}$ of the block matrix with a leading 1.
For real positive definite $S$, the loop hafnian $\operatorname{lhaf}(S,v)$ equals the joint non-centered moment $\mathbb{E}[X_1\cdots X_\mu]$ of Gaussian variables with covariance $S$ and mean $v$.
Both sides of Eq. (4) depend continuously on the entries of $S$ and $v$ near $z=0$, so coefficient extraction remains well-defined even for singular or non-diagonalizable matrices.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

Inference: Because the coefficients are polynomial in $S$ and $v$, Eq. (4) should hold as an identity of formal power series with no convergence or positivity hypotheses; the paper's analytic proof and continuity extension suggest but do not formulate this.
Inference: A numerical check on a non-diagonalizable symmetric matrix such as $S=\begin{pmatrix} i&1\\ 1&-i\end{pmatrix}$ would isolate whether the continuity extension is doing the work; the paper does not perform such a check.
Inference: The derivative machinery connecting loop hafnians to Gaussian exponentials could be reused to derive recurrences for fast loop-hafnian evaluation, though no algorithm is proposed here.

Share X Bluesky LinkedIn Reddit HN

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Desk Editor's Note

read the letter

The one thing you should know: this note proves the loop hafnian generating function of Eq. (4) for all symmetric complex matrices, discarding the quantum-Gaussian-state caveat from Bulmer et al. The proof is sound and the paper does exactly what it claims. The reader's ACCEPT verdict is right, and the stress-test check holds up on reading.

What's actually new: the formula was derived in [1] by quantum-optical methods for covariance matrices with special block symmetries, then used outside that class as a numerically-supported assumption. Tarasov turns that assumption into a theorem. The method is not novel—it adapts the hafnian master theorem template—but the result is new and useful. It removes a caveat in an active line of quantum optical simulation work. It does not change complexity results or open new applications, so the significance is real but contained.

The paper does well: the proof is self-contained, the Gaussian integration steps are standard and correctly executed, and the extension by density/continuity is valid. The algebra in Eqs. (12)-(13) is correct. The lemma is properly cited and the paper gives a useful probabilistic reading in Remark 1. No circularity: the proof does not rely on the restricted earlier derivation. The references are appropriate, including the author's own hafnian master theorem, which is legitimately the template.

Soft spots are minor. The proof starts by assuming S is diagonalizable and invertible with a symmetric square root; the existence of such a root is not proved or cited. It does exist (a polynomial in S can interpolate a square root on the spectrum), so this is a harmless omission. The text says 'any of its square roots, which is also symmetric,' which is slightly careless—not every square root is symmetric—but the proof only needs one, and existence is enough. The analytic continuation step is stated briefly but is standard: both sides are continuous in S and v for small z near the origin, and the derivatives at z=0 pass to the limit. These are wording issues, not mathematical gaps. The odd-dimensional case is handled by a padding trick in Remark 2, which is fine.

Who is this for? Anyone implementing loop hafnian algorithms or GBS simulations who wants a rigorous justification for using the generating function beyond the quantum covariance class. The note is short and clearly written. It deserves a serious referee; the referee's job would be to tighten the square-root wording and the continuity argument, not to doubt the core. I'd bring it to the reading group, and I'd cite it if I were writing on Gaussian boson sampling. Recommend peer review.

Referee Report

0 major / 5 minor

Summary. The manuscript proves that the loop-hafnian generating function identity (Eq. (4)) is valid for every complex symmetric 2m x 2m matrix S and every 2m-dimensional complex vector v. The proof uses the derivative representation of loop hafnians (the Lemma), rewrites the exponential via a Gaussian integral for diagonalizable invertible S (Eq. (9)), converts the x-derivatives into z-derivatives of an auxiliary Gaussian exponent (Eq. (10)), evaluates the Gaussian integral, simplifies the matrix algebra (Eqs. (12)-(13)), and finally extends by continuity to all symmetric S. The note also connects loop hafnians to Gaussian moments (Remark 1) and discusses odd-dimensional matrices (Remark 2).

Significance. If correct, the result removes the quantum-Gaussian-state restrictions of the earlier derivation in [1], making the loop-hafnian master theorem a general identity for arbitrary symmetric complex matrices. The proof is self-contained, elementary, and clearly written; it does not rely on the quantum-optical context and uses only standard Gaussian integration plus a polynomial density argument. This is a useful contribution for Gaussian boson sampling and related combinatorial applications, although the mathematical content is not deep and the final identity is not unexpected given the existing hafnian master theorem.

minor comments (5)

[Proof of the theorem, before Eq. (9)] The sentence 'Let the matrix S^{1/2} be any of its square roots, which is also a diagonalizable and invertible 2m x 2m complex symmetric matrix' is not literally true, since not every square root of a symmetric matrix is symmetric. The argument only needs the existence of a symmetric square root, which follows for diagonalizable invertible S by polynomial interpolation on the spectrum; please rephrase to avoid asserting that every square root has this property.
[Proof of the theorem, final paragraph] The density and continuity extension is compressed into a single sentence. It would be helpful to state explicitly that diagonalizable invertible complex symmetric matrices are dense in the space of all complex symmetric matrices, that for z in a sufficiently small neighborhood of the origin det(1-ZS) is bounded away from zero, and that equality of all z-derivatives at z=0 passes to the limit.
[Theorem, Eq. (4)] The notation \tilde S_{\{n_j\},\{n_j\}} with a doubled subscript is redundant, since the construction depends only on the single tuple \{n_j\}. Using the same notation as in the Lemma (e.g., \tilde S_{\{n_j\}}) would avoid confusion.
[Lemma proof, Eq. (6)] The text introduces the Fa di Bruno calculation as 'a comment following their arguments,' but the paragraph is a complete proof of the lemma. The wording could be changed to 'proof' or 'self-contained derivation' to reflect what is actually presented.
[Theorem, Eq. (4)] The right-hand side of Eq. (4) is an infinite power series. The paper would benefit from a remark that the identity is a formal power series identity, or that for fixed S and sufficiently small |z_j| the series converges and both sides are analytic.

Circularity Check

0 steps flagged · score 0.0 of 10

No significant circularity: the Gaussian-integral proof is self-contained, and the cited prior results are non-load-bearing.

full rationale

The theorem is derived from first principles within the paper. The starting point is the derivative identity in Eq. (5), which is independently justified by the multivariate Faà di Bruno argument in Eqs. (6)-(7) rather than by importing the generating function. The Gaussian representation in Eq. (9) relies on a standard shifted Gaussian integral identity, and the existence of a symmetric square root for a diagonalizable invertible symmetric matrix is a background algebraic fact, not an assumption of the theorem. The extension from diagonalizable invertible matrices to arbitrary symmetric matrices is made by an explicit continuity and density argument, which is a legitimate limiting procedure rather than a circular appeal. Citations to the author's own prior work, especially [7] and [10], serve only as a technical template and historical context; the central matrix algebra, including the simplification in Eq. (13), is carried out in the paper. No fitted parameters appear, no restricted version of the result is assumed, and no uniqueness theorem is invoked to force the conclusion. I therefore find no significant circularity.

Assumptions & free parameters 0 free parameters · 4 assumptions · 0 invented entities

No free parameters are introduced; the proof uses standard Gaussian integration and a density argument. No new mathematical objects are invented; the loop hafnian and its generating function already existed in the literature.

assumptions (4)

standard math Standard Gaussian integral formula for quadratic forms with complex coefficients and positive-definite real part
Used in Eq. (9) and in computing the integral in Eq. (12). The extension to complex S relies on analytic continuation, which is not explicitly discussed.
standard math The set of diagonalizable invertible matrices is dense in the set of complex symmetric matrices, and both sides of Eq. (4) are continuous in S near z=0
Justifies the final extension by continuity from diagonalizable invertible S to all symmetric S. This is a standard density argument.
standard math Every diagonalizable invertible complex symmetric matrix has a diagonalizable invertible symmetric square root
Used to set up the Gaussian representation in Eq. (9). The paper does not cite or prove this, but it is a known result.
standard math Multivariate Faa di Bruno formula for the derivative of a composition
Used in the comment proving the lemma, Eq. (6).

how reviews work

0 comments

Cite this review

Pith. "Pith review of Derivation of the Loop Hafnian Generating Function for Arbitrary Symmetric Matrices via Gaussian Integration." pith.science (2026). https://pith.science/paper/QZCSI6MG

@misc{pith2026250716100,
  author       = {Pith},
  title        = {Pith review of: Derivation of the Loop Hafnian Generating Function for Arbitrary Symmetric Matrices via Gaussian Integration},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/QZCSI6MG}},
  note         = {Machine review of arXiv:2507.16100}
}

read the original abstract

This short note shows that the recently proposed generating function for loop hafnians -- originally derived using quantum-optical methods for a restricted class of matrices -- is in fact valid for arbitrary symmetric matrices. The proof relies solely on Gaussian integration and does not assume any additional properties inherited from the covariance matrices of quantum Gaussian states.

Discussion (0). Sign in to comment.

Reference graph

Works this paper leans on

10 extracted references · 7 canonical work pages

[1]

J. F. Bulmer, J. Mart ´ ınez-Cifuentes, B. A. Bell, and N. Quesada, Simulating lossy and partially distinguishable quan- tum optical circuits: theory, algorithms and applications to experiment validation and state preparation, arXiv preprint: 10.48550/arXiv.2412.17742 (2024)

work page Pith review arXiv doi:10.48550/arxiv.2412.17742 2024
[2]

Bj¨ orklund, B

A. Bj¨ orklund, B. Gupt, and N. Quesada, A faster hafnian formula for complex matrices and its benchmarking on a supercomputer, Journal of Experimental Algorithmics24, 1 (2019)

work page 2019
[3]

Quesada, L

N. Quesada, L. G. Helt, J. Izaac, J. M. Arrazola, R. Shahrokhshahi, C. R. Myers, and K. K. Sabapathy, Simulating realistic non-Gaussian state preparation, Phys. Rev. A100, 022341 (2019)

work page 2019
[4]

Kruse, C

R. Kruse, C. S. Hamilton, L. Sansoni, S. Barkhofen, C. Silberhorn, and I. Jex, Detailed study of Gaussian boson sampling, Phys. Rev. A100, 032326 (2019)

2019
[5]

C. S. Hamilton, R. Kruse, L. Sansoni, S. Barkhofen, C. Silberhorn, and I. Jex, Gaussian Boson Sampling, Phys. Rev. Lett. 119, 170501 (2017)

2017
[6]

Combinatorics of Partial Derivatives

M. Hardy, Combinatorics of partial derivatives, arXiv preprint: 10.48550/arXiv.math/0601149 (2006)

work page Pith review arXiv doi:10.48550/arxiv.math/0601149 2006
[7]

V. V. Kocharovsky, Vl. V. Kocharovsky, and S. V. Tarasov, The Hafnian Master Theorem, Linear Algebra and its Appli- cations651, 144 (2022)

work page 2022
[8]

Isserlis, On a formula for the product-moment coefficient of any order of a normal frequency distribution in any number of variables, Biometrika12, 134 (1918)

L. Isserlis, On a formula for the product-moment coefficient of any order of a normal frequency distribution in any number of variables, Biometrika12, 134 (1918)

1918

Show all 10 references

[9]

Wick, The evaluation of the collision matrix, Phys

G.-C. Wick, The evaluation of the collision matrix, Phys. Rev.80, 268 (1950)

1950
[10]

V. V. Kocharovsky, Vl. V. Kocharovsky, and S. V. Tarasov, Atomic boson sampling in a Bose-Einstein-condensed gas, Phys. Rev. A106, 063312 (2022)

2022

Pith tools

Reviewed August 6, 2026 · model on record in the stance chip above.

[1] [1]

J. F. Bulmer, J. Mart ´ ınez-Cifuentes, B. A. Bell, and N. Quesada, Simulating lossy and partially distinguishable quan- tum optical circuits: theory, algorithms and applications to experiment validation and state preparation, arXiv preprint: 10.48550/arXiv.2412.17742 (2024)

work page Pith review arXiv doi:10.48550/arxiv.2412.17742 2024

[2] [2]

Bj¨ orklund, B

A. Bj¨ orklund, B. Gupt, and N. Quesada, A faster hafnian formula for complex matrices and its benchmarking on a supercomputer, Journal of Experimental Algorithmics24, 1 (2019)

work page 2019

[3] [3]

Quesada, L

N. Quesada, L. G. Helt, J. Izaac, J. M. Arrazola, R. Shahrokhshahi, C. R. Myers, and K. K. Sabapathy, Simulating realistic non-Gaussian state preparation, Phys. Rev. A100, 022341 (2019)

work page 2019

[4] [4]

Kruse, C

R. Kruse, C. S. Hamilton, L. Sansoni, S. Barkhofen, C. Silberhorn, and I. Jex, Detailed study of Gaussian boson sampling, Phys. Rev. A100, 032326 (2019)

2019

[5] [5]

C. S. Hamilton, R. Kruse, L. Sansoni, S. Barkhofen, C. Silberhorn, and I. Jex, Gaussian Boson Sampling, Phys. Rev. Lett. 119, 170501 (2017)

2017

[6] [6]

Combinatorics of Partial Derivatives

M. Hardy, Combinatorics of partial derivatives, arXiv preprint: 10.48550/arXiv.math/0601149 (2006)

work page Pith review arXiv doi:10.48550/arxiv.math/0601149 2006

[7] [7]

V. V. Kocharovsky, Vl. V. Kocharovsky, and S. V. Tarasov, The Hafnian Master Theorem, Linear Algebra and its Appli- cations651, 144 (2022)

work page 2022

[8] [8]

Isserlis, On a formula for the product-moment coefficient of any order of a normal frequency distribution in any number of variables, Biometrika12, 134 (1918)

L. Isserlis, On a formula for the product-moment coefficient of any order of a normal frequency distribution in any number of variables, Biometrika12, 134 (1918)

1918