Gaussian Boson Sampling as a Calibrated Variance-Reduction Engine: A DSFL Explanation

1. Introduction

Gaussian Boson Sampling (GBS) produces photon-count patterns whose probabilities are governed by hafnian functionals of a covariance encoded in a Gaussian photonic state. Beyond its foundational role in quantum advantage, GBS has recently been positioned as a practical engine for approximating Gaussian expectation problems (GEPs), where one seeks $\mu =\mathbb{E}\left[f\right(X\left)\right]$ for $X\sim \mathcal{N}(0,B)$ and analytic or polynomial f. Wick’s theorem expresses $\mu$ as a hafnian–weighted sum over multi-indices drawn from B, aligning the target with the native GBS statistics. Building on this structure, Andersen and Shan introduced degree-resolved GBS estimators and proved open, nonempty regions in which GBS achieves $(\epsilon ,\delta )$ sample–complexity advantages over Monte Carlo (MC), and further showed that tuning the mean photon number to the task increases the percentage of instances with advantage [1,2].

1.1. Calibrated Residuals, Principal–Angle Conditioning, and Predict–Then–Prove Guarantees

This work advances three aspects that are often treated separately. First, the residual $R_K\left(t\right)$ is not merely a distance: it is a calibrated, minimal, and dynamical quantity. We prove that every admissible tuning step is $L^2$–nonexpansive, so the residual cannot inflate under retuning or composition. Second, conditioning is captured quantitatively by the principal angle ${\theta }_F$ (via its cosine $\gamma$), and this geometric control enters directly into the variance and MSE formulae, providing sharp two–sided bounds; this bridges statistical estimation, photonic physics, and numerical analysis in a single comparison geometry. Third, the framework is predict–then–prove: before sampling, it specifies where to operate (photon–number matching), how conditioning will affect dispersion (principal angles), and what decay rate to expect (Lyapunov envelope), and then establishes these prescriptions as theorems.

1.2. Problem

Practitioners need to foresee, before allocating shots, when GBS will outperform MC and how to tune the device to realize that advantage. Concretely: given a covariance B and an integrand f, can we predict in advance whether calibrated GBS will require asymptotically fewer samples than MC to reach a prescribed $(\epsilon ,\delta )$ accuracy, and can we prescribe a robust one-parameter tuning rule that maximizes the gain? We seek design-time criteria (not post hoc curve fits) that tell us when to use GBS, when to defer to MC, and how to switch between them by degree.

1.3. DSFL Intuition: Forecasting When GBS Beats MC

The DSFL framework predicts, a priori, when and how GBS will outperform Monte Carlo (MC) by placing the target statistics and the sampler in a single Hilbert geometry and tracking one quantity—the residual. This yields the following design-time principles:

• Single yardstick. Target and sampler are embedded in the same space; a single residual provides a quantitative knob to optimize and monitor.
• No inflation under tuning. Retuning, counting, and block composition act as nonexpansive maps; the residual cannot increase under admissible operations.
• Unique, certifiable optimum. The dominant variance/MSE factor is strictly convex in the scale and attains its unique minimum at the mean–photon–number match $m_{tB}=2K$; this yields a robust one-line tuning rule.
• Predictable conditioning. Alignment between sampling and blueprint features is quantified by the principal angle ${\theta }_F$ (cosine $\gamma$); small $\gamma$ certifies low dispersion. A computable surrogate from warm-up features allows pre-run assessment.
• Rates, not just limits. If each update removes a fixed fraction of the residual, the residual obeys a non-asymptotic exponential envelope, providing a predicted semi-log MSE slope before data collection.
• Regime map and hybrid policy. The geometry identifies regimes where GBS is inherently favorable (high degrees) and where MC suffices (low degrees), and prescribes a degree-wise hybrid switch using the same residual/angle diagnostics.
• Robustness. Composition and small hardware drifts degrade angles and rates smoothly; re-matching $m_{tB}$ restores near-optimality.

In summary, DSFL turns GBS from a black-box sampler into a calibrated, predictable engine: it specifies where to operate, how conditioning impacts dispersion, and what decay rate to expect, enabling an informed, pre-run decision on when—and by how much—GBS will outperform MC.

1.4. Thesis (DSFL View)

We show that GBS is a calibrated variance-reduction engine inside a single Hilbert comparison geometry governed by a Deterministic Statistical Feedback Law (DSFL). This provides one residual to minimize and four provable pillars: (i) Admissibility & DPI: counting and retuning are $L^2$-nonexpansive and intertwine with the estimator, so the calibrated residual cannot inflate; (ii) Calibration optimality: the unique minimizer of the leading variance/MSE factor is the photon-number match $m_{tB}=2K$, which implements the tuning advocated and empirically optimized in [1,2]; (iii) Angle conditioning: principal (Friedrichs) angles between the sampling and blueprint subspaces yield sharp two-sided dispersion bounds $(1\pm cos{\theta }_F)$, explaining conditioning and leakage; (iv) Lyapunov envelope: a one-step contraction turns into a non-asymptotic exponential error decay with rate $\alpha \left(t\right)$ maximized at the photon match.

1.5. Contributions

• A unified geometry that proves a data–processing inequality (no residual inflation) for tuning and composition, and recovers the $(\epsilon ,\delta )$ envelopes of [1].
• A one-line strict-convexity proof that predicts the optimal tuning rule $m_{tB}=2K$, matching the optimized estimators in [3].
• Principal-angle bounds that forecast dispersion and a Lyapunov inequality that delivers exponential, non-asymptotic decay—together forming practical, pre-run design rules to decide GBS vs. MC and to set t.
• A regime map that explains, a priori, when calibrated GBS is polynomial while MC is exponential (large degree), and a hybrid, degree-wise policy for the small-degree regime.

In short, DSFL turns “folk rules’’ into theorems: match $m_{tB}$, measure angles, and expect exponential envelopes—thereby foreseeing and tuning the regimes where GBS outperforms MC, in line with and extending the results of [1,2].

2. What GBS Does (Succinct but Precise)

GBS distribution and hafnians 

Let $B\in {\mathbb{R}}^{N\times N}$ be symmetric with $0<\mathrm{spec}\left(B\right)<1$. The GBS probability of a pattern $I\in {\mathbb{N}}^N$ is

$$P_B\left(I\right)={\displaystyle \frac{d_B}{I!}}Haf{\left(B_I\right)}^2,d_B:=\prod _{j=1}^N\sqrt{1-{\lambda }_j^2},$$

(1)

where $\left\{{\lambda }_j\right\}$ are eigenvalues of B. The mean photon number is

$$m_B=\sum _{j=1}^N{\displaystyle \frac{{\lambda }_j^2}{1-{\lambda }_j^2}}.$$

(2)

Gaussian expectations via Wick’s theorem 

For $f\left(x\right)={\sum }_{k\ge 0}{\sum }_{\left|I\right|=k}{a}_I{x}^I$ and $X\sim \mathcal{N}(0,B)$,

$$\mu =\mathbb{E}\left[f\left(X\right)\right]=\sum _{k\ge 0}\sum _{\left|I\right|=2k}{a}_{I}Haf\left(B_I\right).$$

(3)

A degree partition gives $\mu ={\sum }_{k=0}^K{\mu }_k$ with ${\mu }_k={\sum }_{\left|I\right|=2k}{a}_{I}Haf\left(B_I\right)$.

Two practical GBS estimators 

GBS-I (importance). For ${\mu }_k^{\left(2\right)}={\sum }_{\left|I\right|=2k}{a}_{I}Haf{\left(B_I\right)}^2$, the estimator using weights $I!a_I/d_B$ is unbiased.

GBS-P (probability). For ${\mu }_k={\sum }_{\left|I\right|=2k}{a}_{I}Haf\left(B_I\right)$, estimate $p_J=P_B\left(J\right)$ by counts and average $\sqrt{p_J}$ with weights. The bias vanishes at rate $1/n$; the leading term is explicit.

Tuning knob: photon-number matching 

Scale $B\mapsto tB$ with $t\in (0,{\lambda }_1^{-1})$ and choose t so that

$$m_{tB}=\sum _{j=1}^N{\displaystyle \frac{t^2{\lambda }_j^2}{1-t^2{\lambda }_j^2}}=2K,$$

(4)

i.e. match the sampler’s mean photon number to the target degree K.

3. Related Work

This section situates our DSFL-based treatment of Gaussian Boson Sampling (GBS) within three intertwined literatures: (i) the foundations of GBS, including sampling distributions, computational hardness, and experimental progress; (ii) GBS as a vehicle for numerical estimation, in particular its use for Gaussian expectation problems together with formal $(\epsilon ,\delta )$-sample–complexity claims and percentage-of-advantage quantification; and (iii) control-theoretic and geometric principles that justify our residual-contraction, angle-conditioned bounds, and data–processing inequalities in a single Hilbertian comparison geometry. We use these threads to explain how the “folk rules’’ of GBS operation—photon-number tuning, conditioning via feature alignment, and non–asymptotic error envelopes—can be elevated to theorems, and to place our results alongside prior theory and practice.

3.1. Foundations of GBS and Hardness

The modern formulation of GBS as a sampling task driven by hafnian-valued amplitudes is due to Hamilton et al. [4] with detailed modeling and refinements in [5]. Complexity-theoretic hardness follows the linear-optical lineage of Aaronson–Arkhipov [6] and extends to experimentally informed, high-dimensional variants [7,8]. Algorithmic primitives for hafnian evaluation and classical simulation—including supercomputer-grade implementations and the thewalrus library—are documented in [9,10]. Photonic demonstrations approaching or attaining quantum advantage are reported in [11,12]. Related Gaussian models include loop-hafnians for displaced states [13] and threshold-detection (Torontoian) statistics [14], which frame hardware-relevant variants. On the application side, GBS has been explored for vibronic spectra [15,16], molecular docking [17], dense subgraph discovery and graph kernels [18,19,20,21], random point processes [22,23], and scalable time-bin/loop architectures [24,25]. These works supply the physical and computational substrate on which our estimation-theoretic analysis is built.

3.2. GBS Estimators for Gaussian Expectations

The use of GBS samples to approximate Gaussian expectation problems (GEPs)—via Wick’s theorem [26] to express Gaussian moments as hafnian-weighted functionals—was developed in [1], which introduced two concrete estimators (GBS-I and GBS-P) and proved non-asymptotic $(\epsilon ,\delta )$ guarantees on an open, nonempty subset of instances. An optimized variant matching the device’s mean photon number to the target degree and quantifying the percentage of the problem space with GBS advantage appears in [2]. Our DSFL framework recovers these bounds within a single comparison geometry, explains the tuning law ($m_{tB}=2K$), and adds angle-conditioned rate controls. Complementary analytic identities (e.g., the hafnian master theorem [27]) and approximation tools [28,29] delineate tractable regimes and provide context for the sampling-based estimators considered here.

3.3. Hilbertian Residual Control and Data Processing

Our analysis leverages standard Hilbert-space contractivity of orthogonal projections (counting) and parameter changes (retuning) to derive a data–processing inequality (DPI) for the DSFL residual, together with principal-angle (Friedrichs) bounds for conditioning and Lyapunov-type decay envelopes. While these tools are classical (see, e.g., two-subspace/Friedrichs-angle inequalities and discrete Grönwall arguments), their joint instantiation for GBS estimation and calibration appears to be new in this context. For Monte Carlo baselines and variance-reduction perspectives that we compare against, see [30,31,32,33].

Gaussian Boson Sampling (GBS) was introduced as a photonic sampling model whose output probabilities are governed by hafnians of submatrices associated with a Gaussian covariance [4,5]. In analogy with (standard) Boson Sampling [6], a growing body of work establishes classical intractability (under plausible conjectures) not only for exact sampling but also for sampling from suitable approximations to the ideal distribution [7,8]. Algorithmic progress on hafnian evaluation—both general-purpose [9] and software implementations (e.g. The Walrus) [10]—has enabled faithful classical simulations at nontrivial scales, supporting experimental validation [11,12] and methodological development.

3.4. Applications of GBS

Proposed applications of GBS span vibronic spectroscopy [15,16], molecular docking [17], dense subgraph discovery [18], perfect matchings [19], graph similarity and kernels [20,21], point processes [22], and photonic simulations of molecular dynamics [23]. These lines have been complemented by platform-level studies of time-bin/loop architectures that aim at scalable implementations [24,25].

3.5. GBS Estimators for Gaussian Expectations

The specific connection between GBS and Gaussian expectation problems (GEPs) is made precise by Wick’s theorem [26], which expresses Gaussian moments as hafnians. Building on this structure, Andersen and Shan (i) introduced degree-resolved GBS estimators (GBS-P and GBS-I) for GEPs and proved open, nonempty regions of exponential advantage over plain Monte Carlo (MC), measured in guaranteed $(\epsilon ,\delta )$ sample sizes [1], and (ii) optimized the mean photon number to the task and quantified the fraction of the problem space displaying advantage, with large observed percentages and sizeable efficiency ratios [2]. Related algorithmic components include hafnian bounds (via Stirling-type inequalities) and polylogarithmic envelopes that deliver computable, nonasymptotic sample-complexity expressions.

3.6. Variance-Reduction and MC Baselines

From a classical perspective, GBS-I can be seen as an importance sampler whose proposal is the physical GBS law; GBS-P corresponds to a (biased-but-consistent) probability estimator with explicit leading error terms. These sit alongside standard MC variance-reduction techniques (importance, stratification, antithetic sampling, quasi-MC), whose strengths and limitations in high dimensions are well documented [30,31,32,33]. The present analysis adopts plain MC as the baseline, following [1,2].

3.7. Hafnian Identities and Tractable Limits

Structural identities such as the hafnian master theorem [27] provide tractable regimes for particular coefficient ensembles (e.g. $z^I/I!$ weights), while general approximation results for hafnians and permanents (e.g. Barvinok-type bounds) furnish algorithmic surrogates outside the photonic setting [28,29]. These are complementary to the sampling-based estimators considered here.

3.8. Angles, Data Processing, and Residual Control

The present DSFL view organizes the GBS estimation task in one Hilbert geometry and leverages two ingredients that are ubiquitous in information processing and operator theory: (i) data–processing inequalities (DPI) for nonexpansive, intertwining updates (counting maps, retuning), ensuring no residual inflation under admissible operations; and (ii) principal-angle (Friedrichs) bounds that control conditioning of two intertwined subspaces (blueprint vs sampling features) and yield two-sided variance inequalities with factors $(1\pm cos{\theta }_F)$. The Lyapunov envelope for the residual—a discrete coercivity statement converting one-step contraction into an exponential decay rate—is classical in stability theory and is instantiated here by the t-matched (photon-number–matched) tuning law. While the DSFL terminology is specific to this work, the underlying tools are aligned with standard Hilbertian projection theory and CS-decompositions.

3.9. Positioning

Relative to [3,3], the present manuscript provides a unifying explanation of when/why/how GBS outperforms MC by: (a) proving that mean–photon–number matching is the unique minimizer of the leading residual/variance factor; (b) quantifying conditioning via principal angles with explicit two-sided bounds; and (c) deriving nonasymptotic Lyapunov envelopes for error decay that recover the $(\epsilon ,\delta )$ sample-size laws. This framework is compatible with current hardware models [34,35,36] and classical tooling for hafnians [9,10], and it complements domain applications [15,17,18,19,20,22,23].

4. Limitations and Scope

The guarantees developed here are tailored to Gaussian expectation problems in which Wick reduction expresses the target as hafnian–weighted functionals of a real, symmetric covariance with spectrum strictly inside $(-1,1)$. This setting aligns the photonic sampler’s distribution with the blueprint geometry in a way that makes DSFL calibration, principal-angle conditioning, and Lyapunov envelopes directly applicable. Extending the framework beyond this regime entails additional structure. For non-Gaussian observables or general boson sampling models (e.g., loop-hafnians with coherent displacements, threshold detection/Torontoian statistics, or complex-valued amplitudes with sign patterns), one must re-specify the feature map, the admissible updates, and the calibration functional so that (i) an appropriate inner-product geometry captures the relevant moments, (ii) the data–processing inequality remains valid under the physical update semigroup, and (iii) bias terms induced by nonlinearity or signs are explicitly controlled. A second limitation concerns conditioning proxies: our use of principal-angle surrogates (e.g., cross-Gram norms between sampling and blueprint features) is computationally light but approximate; computing exact projectors or Friedrichs angles may be prohibitive for very large $(N,K)$. Finally, while our analysis covers noiseless sampling and block composition via nonexpansive maps, hardware noise and calibration drift require model-specific contractions; the DSFL conclusions persist under such perturbations only insofar as a positive coercivity margin can be certified after re-matching the mean photon number.

4.1. DSFL Foresight for GBS: Verbal Logic → Mathematical Predictors

We give a from-first-principles account of how the Deterministic Statistical Feedback Law (DSFL) predicts (before any run) the behavior of Gaussian Boson Sampling (GBS): how tuning affects error, when/why it beats Monte Carlo (MC), how variance/bias scales, why mean-photon-number matching is optimal, and how composition/noise impact performance. Each verbal claim is paired with an explicit, checkable mathematical statement.

5. Stage 1: Put GBS in One Hilbert Geometry (the DSFL Stage)

DSFL places the statistical blueprint (the quantity to estimate) and the physical sampler (the device) in the same comparison geometry and measures a single misfit (residual). In our setting, GBS turns a scale parameter t (the photon-number knob) into a sample law $P_{tB}$. The estimator (GBS-I or GBS-P restricted to degree K) is a calibration map ${\mathcal{I}}_t$ that reads off ${\mu }_K$ from $P_{tB}$. Alignment ⇒ small residual; misalignment ⇒ large residual.

Mathematics:

Fix a degree K. Let the blueprint vector be

$$s_K\equiv {\left(a_I\right)}_{\left|I\right|=2K}\in {\mathbb{R}}^{D_K},D_K=\#\{I:\left|I\right|=2K\}.$$

(5)

Let the physical law at scale t be the degree-$2K$ marginal

$$p^{\left(t\right)}={\left(p_I^{\left(t\right)}\right)}_{\left|I\right|=2K},p_I^{\left(t\right)}={\displaystyle \frac{d_{tB}}{I!}}Haf{\left({\left(tB\right)}_I\right)}^2,d_{tB}=\prod _{j=1}^N\sqrt{1-t^2{\lambda }_j^2}.$$

(6)

Work in $\mathcal{H}={\mathbb{R}}^{D_K}$ with the standard inner product. Define the readout (calibration) ${\mathcal{I}}_t:\mathcal{H}\to \mathcal{H}$ by the chosen estimator:

$$\begin{array}{cc}\hfill \text{GBS-P}\left(\mathrm{probability}\right)\text{:}& {\left({\mathcal{I}}_{t}p\right)}_K\approx \sum _{\left|I\right|=2K}{a}_I{\displaystyle \frac{p_I^{\left(t\right)}}{d_{tB}}},\hfill \\ \hfill \text{GBS-I}(\mathrm{importance}\mathrm{on}{Haf}^2)\text{:}& {\left({\mathcal{I}}_{t}p\right)}_K\approx \sum _{\left|I\right|=2K}{a}_I{\displaystyle \frac{I!}{d_{tB}}}{p}_I^{\left(t\right)}.\hfill \end{array}$$

Define the Residual of Sameness

$${\mathcal{R}}_K\left(t\right)\equiv \parallel p^{\left(t\right)}-{\mathcal{I}}_t{s}_K{\parallel }_{\mathcal{H}}^2.$$

(7)

This single scalar is the DSFL misfit we will minimize, bound, and track.

6. DSFL Primer (Minimal Machinery, with Proofs)

6.1. Comparison Geometry

Let $(\mathcal{H},\langle ·,·\rangle )$ be a (real or complex) Hilbert space equipped with the norm $\parallel x\parallel :=\sqrt{\langle x,x\rangle }$. We fix two closed subspaces $\mathsf{S},\mathsf{P}\subset \mathcal{H}$ (“statistical” and “physical” channels), and denote by $P_{\mathsf{S}},P_{\mathsf{P}}$ the corresponding orthogonal projectors.

6.2. Calibration (Interchangeability)

Definition 6.1 

(Calibration pair). A calibration is a pair of bounded operators

$$\mathcal{I}:\mathsf{S}\to \mathsf{P},\mathcal{J}:\mathsf{P}\to \mathsf{S},$$

(29)

satisfying the interchangeability identities

$$\mathcal{I}\mathcal{J}={id}_{\mathsf{P}},\mathcal{J}\mathcal{I}=P_{\mathsf{S}}.$$

(30)

Remark 6.1. 

The relations (30) imply that $\mathcal{I}$ is a left inverse of $\mathcal{J}$ on $\mathsf{P}$, and $\mathcal{J}$ is a right inverse of $\mathcal{I}$ on $\mathsf{S}$. In particular, $\mathcal{I}$ is an isomorphism of $\mathsf{S}/ker\mathcal{I}$ onto $\mathsf{P}$, and $\mathcal{J}$ restricts to the identity on $\mathcal{I}\left(\mathsf{S}\right)$.

6.3. Residual of Sameness

Definition 6.2 

(Residual). Given $s\in \mathsf{S}$ (blueprint) and $p\in \mathsf{P}$ (response), the Residual of Sameness is the squared $\mathcal{H}$–distance

$$\mathcal{R}(s,p):={\parallel p-\mathcal{I}s\parallel }_{\mathcal{H}}^2.$$

(31)

Remark 6.2. 

By construction, $\mathcal{R}(s,p)=0$ iff $p=\mathcal{I}s$. Moreover, $\mathcal{R}$ is invariant under any isometry U on $\mathcal{H}$ that preserves $\mathrm{ran}\left(\mathcal{I}\right)$: if $U\left(\mathrm{ran}\mathcal{I}\right)\subseteq \mathrm{ran}\mathcal{I}$, then $\parallel Up-U\mathcal{I}s\parallel =\parallel U(p-\mathcal{I}s)\parallel =\parallel p-\mathcal{I}s\parallel$.

6.4. Admissibility and the Data–Processing Inequality

Definition 6.3 

(Admissible update). A pair of bounded operators $(\tilde{\mathrm{\Phi }},\mathrm{\Phi })$ with

$$\tilde{\mathrm{\Phi }}:\mathsf{S}\to \mathsf{S},\mathrm{\Phi }:\mathsf{P}\to \mathsf{P},$$

(32)

is admissible if:

1.

Intertwining (calibration consistency):

$$\mathrm{\Phi }\mathcal{I}=\mathcal{I}\tilde{\mathrm{\Phi }};$$

(33)

2.

Nonexpansiveness: ${\parallel \mathrm{\Phi }\parallel }_{\mathcal{H}\to \mathcal{H}}\le 1$.

Theorem 6.4 

(Data–processing inequality (DPI)). If $(\tilde{\mathrm{\Phi }},\mathrm{\Phi })$ is admissible, then for all $s\in \mathsf{S}$ and $p\in \mathsf{P}$,

$$\mathcal{R}(\tilde{\mathrm{\Phi }}s,\mathrm{\Phi }p)\le \mathcal{R}(s,p).$$

(34)

Proof. 

Using (33),

$$\mathcal{R}(\tilde{\mathrm{\Phi }}s,\mathrm{\Phi }p)=\parallel \mathrm{\Phi }p-\mathcal{I}\tilde{\mathrm{\Phi }}s{\parallel }^2={\parallel \mathrm{\Phi }p-\mathrm{\Phi }\mathcal{I}s\parallel }^2={\parallel \mathrm{\Phi }(p-\mathcal{I}s)\parallel }^2\le {\parallel \mathrm{\Phi }\parallel }^2\parallel {p-\mathcal{I}s\parallel }^2\le {\parallel p-\mathcal{I}s\parallel }^2,$$

(35)

which is (34).    □

Corollary 6.5 

(Iterated DPI). If ${\left\{({\tilde{\mathrm{\Phi }}}_k,{\mathrm{\Phi }}_k)\right\}}_{k=1}^m$ are admissible, then so is $(\tilde{\mathrm{\Phi }},\mathrm{\Phi })=({\prod }_{k=1}^m{\tilde{\mathrm{\Phi }}}_k,{\prod }_{k=1}^m{\mathrm{\Phi }}_k)$, and

$$\mathcal{R}\left(\tilde{\mathrm{\Phi }}s,\mathrm{\Phi }p\right)\le \mathcal{R}(s,p).$$

(36)

Proof. 

Intertwining and $\parallel ·\parallel \le 1$ are preserved under composition. Apply Theorem 6.4.    □

6.5. Angles and Conditioning

Let $M,N\subset \mathcal{H}$ be closed subspaces with orthogonal projectors $P_M,P_N$.

Definition 6.6 

(Friedrichs angle). The Friedrichs cosine is $\gamma (M,N):=\parallel P_M{P}_N\parallel$, and the Friedrichs angle ${\theta }_F(M,N)\in [0,\frac{\pi }{2}]$ is defined by

$$cos{\theta }_F(M,N):=\gamma (M,N).$$

(37)

We recall a standard two–subspace inequality; we include a proof for completeness.

Proposition 6.7 

(Two–subspace bounds). Let $M,N\subset \mathcal{H}$ be closed subspaces and set $\gamma =\gamma (M,N)$. For any $a\in M$, $b\in N$,

$$(1-\gamma )\left({\parallel a\parallel }^2+{\parallel b\parallel }^2\right)\le {\parallel a+b\parallel }^2\le (1+\gamma )\left({\parallel a\parallel }^2+{\parallel b\parallel }^2\right).$$

(38)

Proof. 

By polarization,

$${\parallel a+b\parallel }^2={\parallel a\parallel }^2+{\parallel b\parallel }^2+2\mathrm{Re}\langle a,b\rangle .$$

(39)

Since $a\in M$, $b\in N$, we have $\left|\langle a,b\rangle \right|=|\langle a,P_{M}b\rangle |\le \parallel a\parallel \parallel P_{M}b\parallel \le \gamma \parallel a\parallel \parallel b\parallel$. Thus

$$-2\gamma \parallel a\parallel \parallel b\parallel \le 2\mathrm{Re}\langle a,b\rangle \le 2\gamma \parallel a\parallel \parallel b\parallel .$$

(40)

Using $2\parallel a\parallel \parallel b\parallel \le {\parallel a\parallel }^2+{\parallel b\parallel }^2$ and $-2\parallel a\parallel \parallel b\parallel \ge -(\parallel a{\parallel }^2+\parallel b{\parallel }^2)$ yields (38).    □

Proposition 6.8 

(Through–subspace contraction). Let $M,N\subset \mathcal{H}$ be closed subspaces with $\gamma =\parallel P_M{P}_N\parallel$. Then for any bounded $T:\mathcal{H}\to \mathcal{H}$,

$$\parallel P_{M}T{P}_N\parallel \le \parallel T\parallel \gamma .$$

(41)

In particular, if $\parallel T\parallel \le 1$, then $\parallel P_{M}T{P}_N\parallel \le \gamma =cos{\theta }_F(M,N)$.

Proof. 

By the dual characterization of the operator norm,

$$\parallel P_{M}T{P}_N\parallel =\underset{\parallel x\parallel =\parallel y\parallel =1}{sup}|\langle P_{M}T{P}_{N}x,y\rangle |=\underset{\parallel x\parallel =\parallel y\parallel =1}{sup}|\langle T{P}_{N}x,P_{M}y\rangle |\le \parallel T\parallel \underset{\parallel x\parallel =\parallel y\parallel =1}{sup}\left|\langle P_{N}x,P_{M}y\rangle \right|.$$

(42)

But ${sup}_{\parallel x\parallel =\parallel y\parallel =1}|\langle P_{N}x,P_{M}y\rangle |=\parallel P_M{P}_N\parallel$ (take the singular value decomposition of $P_M{P}_N$). Hence (41).    □

Remark 6.3 

(Interpretation). In applications, M and N are the “pointer” subspaces representing, respectively, the sampling polarization and the blueprint (degree–K) polarization. The factor $\gamma =cos{\theta }_F$ controls leakage and conditioning across the seam: estimates that pass through an $L^2$–contractive kernel are suppressed by at least $cos{\theta }_F$.

6.6. Lyapunov Envelopes (Continuous and Discrete)

We formalize a standard energy–dissipation mechanism that will underlie non–asymptotic error envelopes.

Lemma 6.9 

(Differential energy inequality). Let $e\left(t\right)\in \mathcal{H}$ be (weakly) differentiable and set $R\left(t\right):={\parallel e\left(t\right)\parallel }^2$. Assume there exists a bounded self–adjoint $K=K^{*}⪰0$ and a remainder $g\left(t\right)\in \mathcal{H}$ such that

$${\displaystyle \frac{d}{dt}}R\left(t\right)=-2\langle Ke\left(t\right),e\left(t\right)\rangle +2\langle e\left(t\right),g\left(t\right)\rangle .$$

(43)

If, for some $\kappa >0$ and $\epsilon \ge 0$,

$$\langle Ke,e\rangle \ge {\kappa \parallel e\parallel }^2,\left|\langle e,g\rangle \right|\le {\epsilon \parallel e\parallel }^2\mathit{for}\mathit{all}e\in \mathcal{H},$$

(44)

then

$${\displaystyle \frac{d}{dt}}R\left(t\right)\le -2(\kappa -\epsilon )R\left(t\right).$$

(45)

Proof. 

Combine (43) and (44):

$${\displaystyle \frac{d}{dt}}R\left(t\right)\le -{2\kappa \parallel e\left(t\right)\parallel }^2+2\epsilon {\parallel e\left(t\right)\parallel }^2=-2(\kappa -\epsilon )R\left(t\right).$$

(46)

   □

Theorem 6.10 

(Continuous Lyapunov envelope). Under the hypotheses of Lemma 6.9 with $\alpha :=2(\kappa -\epsilon )>0$,

$$R\left(t\right)\le e^{-\alpha (t-t_0)}R\left(t_0\right)\mathit{for}\mathit{all}t\ge t_0.$$

(47)

Proof. 

Apply Grönwall’s inequality to (45).    □

A discrete analogue covers sampling iterations.

Proposition 6.11 

(Discrete Lyapunov envelope). Let ${\left(R_n\right)}_{n\in \mathbb{N}}$ be nonnegative numbers such that for some $\alpha \in (0,1]$,

$$R_{n+1}\le (1-\alpha )R_n\mathit{for}\mathit{all}n\in \mathbb{N}.$$

(48)

Then

$$R_n\le {(1-\alpha )}^n{R}_0\le e^{-\alpha n}{R}_0.$$

(49)

Proof. 

Iterate (48) and use ${(1-\alpha )}^n\le e^{-\alpha n}$.    □

Remark 6.4 

(Connection to admissibility). In many applications (e.g. residuals computed after one sampling/update step), admissibility of the update and a sectorial coercivity provide (48) in expectation, yielding an expected exponential envelope of the form $\mathbb{E}\left(R_n\right)\le e^{-\alpha n}\mathbb{E}\left(R_0\right)$.

7. Mapping GBS ⟺ DSFL

Here we place the GBS estimation problem and the DSFL control law into a single, explicit Hilbert–space model and show that the main “folk rules” for operating a Gaussian Boson Sampler—how to tune the mean photon number, when advantage over Monte Carlo (MC) appears, how conditioning depends on feature alignment, and why errors decay at a fixed rate—follow as theorems from one comparison geometry. The construction has three ingredients. First, for a fixed target degree K, we embed both the statistical blueprint (the degree–K coefficient vector of the integrand) and the physical sampling law induced by the device at scale t into the same finite–dimensional Hilbert space ${\mathcal{H}}_K={\mathbb{R}}^{D_K}$ indexed by multi–indices ${\mathcal{I}}_K=\{I:|I|=2K\}$. This makes the GBS law $p\left(t\right)$ and the blueprint $s_K$ comparable as vectors and allows us to measure a single misfit—the DSFL residual

$${\mathcal{R}}_{\mathrm{res}}^{\left(K\right)}\left(t\right)={∥p\left(t\right)-\mathrm{\Pi }{s}_K∥}_{{\mathcal{H}}_K}^2.$$

(50)

Second, we model the estimators (GBS–P and GBS–I restricted to degree K) as bounded linear functionals $\mathcal{I}:{\mathcal{H}}_K\to \mathbb{R}$, so that “reading off” the degree–K contribution to the Gaussian expectation is nothing but a linear evaluation in ${\mathcal{H}}_K$. In this form, the basic DSFL admissibility hypotheses hold for the two primitive operations available on hardware or in simulation: (a) changing the scale t (retuning the source) and (b) replacing the law by its empirical counterpart (counting). Both operations act by contractions on ${\mathcal{H}}_K$ and intertwine with the calibration $\mathcal{I}$. Consequently, the DSFL data–processing inequality applies directly: no admissible update can inflate the degree–K residual.

Third, we extract quantitative laws from the geometry. A short calculus identity shows that the leading t–dependence of the residual and of the estimator variance is governed by the convex factor

$$C_{2K}\left(t\right)={\displaystyle \frac{t^{-2K}}{d_{tB}}},\mathrm{with}{\displaystyle \frac{d}{dt}}log{C}_{2K}\left(t\right)={\displaystyle \frac{m_{tB}-2K}{t}},$$

(51)

so the unique minimizer occurs at the photon–number match $m_{tB}=2K$. This is the precise DSFL formulation of “mean–photon–number matching is optimal.” Conditioning and leakage across degrees are controlled by the Friedrichs principal angle ${\theta }_F\left(t\right)$ between the blueprint subspace and the sampling feature subspace in ${\mathcal{H}}_K$; standard two–subspace estimates yield two–sided variance (or squared–bias) bounds with multiplicative factors $(1\pm cos{\theta }_F\left(t\right))$. Finally, a one–step conditional contraction of the residual—guaranteed by admissibility plus a positive angle margin—produces a Lyapunov inequality and hence a non–asymptotic exponential envelope

$$\mathbb{E}{\mathcal{R}}_{\mathrm{res}}^{\left(K\right)}\left(p^{\left(n\right)}\left(t\right)\right)\le e^{-\alpha \left(t\right)n}{\mathcal{R}}_{\mathrm{res}}^{\left(K\right)}\left(p^{\left(0\right)}\left(t\right)\right),$$

(52)

with $\alpha \left(t\right)$ maximal at the photon–match.

The subsections below make this scheme explicit. We first specify the blueprint $s_K$, the physical channel $p\left(t\right)$, and the linear estimator $\mathcal{I}$ inside ${\mathcal{H}}_K$. We then verify admissibility (intertwining and nonexpansiveness) for counting and retuning, yielding a degree–K data–processing inequality. Next, we prove calibration optimality—photon–number matching minimizes the residual—and derive principal–angle bounds for the estimator’s dispersion. We conclude with a Lyapunov–type result that converts the one–step contraction into an exponential decay rate, thereby recovering and clarifying the $(\epsilon ,\delta )$ sample–complexity laws for GBS and the observed advantage regimes over MC.

7.1. Blueprint, Physical Channel, and Calibration in One Geometry

Fix an integer $K\ge 0$ and a symmetric covariance $B\in {\mathbb{R}}^{N\times N}$ with spectrum $0<{\lambda }_N\le \cdots \le {\lambda }_1<1$. Let

$${\mathcal{I}}_K:=\{I=(i_1,\dots ,i_N)\in {\mathbb{N}}^N:\left|I\right|:=i_1+\cdots +i_N=2K\},D_K:=\#{\mathcal{I}}_K.$$

(53)

7.1.1. Hilbert Model

We work on the comparison Hilbert space

$${\mathcal{H}}_K:={\mathbb{R}}^{D_K},\langle x,y\rangle :=\sum _{I\in {\mathcal{I}}_K}{x}_I{y}_I,{\parallel x\parallel }^2:=\sum _I{x}_I^2.$$

(54)

7.1.2. Blueprint

Let $f\left(x\right)={\sum }_{k\ge 0}{\sum }_{\left|I\right|=k}{a}_I{x}_I$ be the target integrand. The degree-K blueprint is the coefficient vector

$$s_K:={\left(a_I\right)}_{\left|I\right|=2K}\in {\mathcal{H}}_K.$$

(55)

Wick’s theorem yields the degree-K contribution to the Gaussian expectation

$${\mu }_K=\sum _{\left|I\right|=2K}{a}_{I}Haf\left(B_I\right),{\mu }_K^{\left(2\right)}=\sum _{\left|I\right|=2K}{a}_{I}Haf{\left(B_I\right)}^2,$$

(56)

where $B_I$ denotes the (multi-)submatrix of B indexed by I.

7.1.3. Physical Channel

For any scale $t\in (0,{\lambda }_1^{-1})$, the GBS distribution on ${\mathcal{I}}_K$ is

$$p_I\left(t\right)={\displaystyle \frac{d_{tB}}{I!}}Haf{\left({\left(tB\right)}_I\right)}^2,d_{tB}=\prod _{j=1}^N\sqrt{1-t^2{\lambda }_j^2}.$$

(57)

We regard $p\left(t\right):={\left(p_I\left(t\right)\right)}_{\left|I\right|=2K}\in {\mathcal{H}}_K$ as the “physical state” at scale t. The (global) mean photon number of $tB$ is

$$m_{tB}:=\sum _{j=1}^N{\displaystyle \frac{t^2{\lambda }_j^2}{1-t^2{\lambda }_j^2}}.$$

(58)

7.1.4. Calibration (Estimator Readout)

We model the degree-K estimators as linear functionals on ${\mathcal{H}}_K$:

$$\begin{array}{cc}\hfill \text{GBS-P}\left(\mathrm{probability}\right)\text{:}& {\mathcal{I}}^{\mathrm{P}}\left(x\right):=\sum _{\left|I\right|=2K}{a}_I{\displaystyle \frac{x_I}{d_{tB}}}=:\langle w^{\mathrm{P}},x\rangle ,w_I^{\mathrm{P}}={\displaystyle \frac{a_I}{d_{tB}}},\hfill \\ \hfill \text{GBS-I}\left(\mathrm{importance}\right)\text{:}& {\mathcal{I}}^{\mathrm{I}}\left(x\right):=\sum _{\left|I\right|=2K}{a}_I{\displaystyle \frac{I!}{d_{tB}}}{x}_I=:\langle w^{\mathrm{I}},x\rangle ,w_I^{\mathrm{I}}={\displaystyle \frac{I!a_I}{d_{tB}}}.\hfill \end{array}$$

Thus $\mathcal{I}:{\mathcal{H}}_K\to \mathbb{R}$ is bounded (by Cauchy–Schwarz). For empirical laws $\widehat{p}\left(t\right)$ obtained by counting, we use the same $\mathcal{I}$; admissibility will follow from linearity and nonexpansiveness of orthogonal projections.

7.1.5. Residual

With this embedding, the DSFL residual for the degree-K task at scale t is

$${\mathcal{R}}_{\mathrm{res}}^{\left(K\right)}\left(t\right):={∥p\left(t\right)-\mathrm{\Pi }{s}_K∥}_{{\mathcal{H}}_K}^2.$$

(59)

where $\mathrm{\Pi }:\mathsf{S}\to {\mathcal{H}}_K$ is the canonical identification $\mathrm{\Pi }{s}_K=s_K$ (we suppress the statistics subspace and view $s_K$ as an element of ${\mathcal{H}}_K$).

7.2. Admissibility of GBS Tuning and Counting

Lemma 7.1 

(Intertwining). Let $\tilde{\mathrm{\Phi }}={id}_{{\mathcal{H}}_K}$ and Φ be either: (i) the identity (fixed t) followed by the orthogonal projection $\mathrm{Proj}$ onto the simplex of probability vectors (counting), or (ii) the change of scale $t\mapsto t^{\prime }$ followed by $\mathrm{Proj}$. Then for $\mathcal{I}\in \{{\mathcal{I}}^{\mathrm{P}},{\mathcal{I}}^{\mathrm{I}}\}$ one has

$$\mathcal{I}\circ \tilde{\mathrm{\Phi }}=\mathcal{I},\mathrm{\Phi }\circ \mathrm{\Pi }=\mathrm{\Pi },\mathit{hence}\mathrm{\Phi }\mathcal{I}=\mathcal{I}\tilde{\mathrm{\Phi }}.$$

(60)

Proof. 

$\tilde{\mathrm{\Phi }}=id$ on the blueprint space; $\mathrm{\Phi }$ acts on physical vectors only and preserves the ambient space ${\mathcal{H}}_K$ (projection and change of parameters do not move $s_K$). Therefore $\mathrm{\Phi }\circ \mathrm{\Pi }=\mathrm{\Pi }$ and $\mathcal{I}\circ \tilde{\mathrm{\Phi }}=\mathcal{I}$, yielding the intertwining $\mathrm{\Phi }\mathcal{I}=\mathcal{I}\tilde{\mathrm{\Phi }}$.    □

Lemma 7.2 

(Nonexpansiveness). The operators Φ of the previous lemma satisfy ${\parallel \mathrm{\Phi }\parallel }_{{\mathcal{H}}_K\to {\mathcal{H}}_K}\le 1$.

Proof. 

Orthogonal projection onto a closed convex set is a contraction in ${\mathcal{H}}_K$; parameter change $t\mapsto t^{\prime }$ is the identity on the ambient linear space (we reparametrize the same coordinates). Hence the composition has operator norm $\le 1$.    □

Proposition 7.3 

(DPI for GBS tuning and counting). Let $(\tilde{\mathrm{\Phi }},\mathrm{\Phi })$ be as above. Then the DPI holds for the degree-K residual:

$${\mathcal{R}}_{\mathrm{res}}^{\left(K\right)}\left(\tilde{\mathrm{\Phi }}{s}_K,\mathrm{\Phi }p\right)\le {\mathcal{R}}_{\mathrm{res}}^{\left(K\right)}(s_K,p)\mathit{for}\mathit{all}{s}_K,p\in {\mathcal{H}}_K.$$

(61)

Proof. 

Apply the abstract DPI to $(\tilde{\mathrm{\Phi }},\mathrm{\Phi })$, with $\mathcal{I}$ as specified and $\mathcal{H}={\mathcal{H}}_K$.    □

7.3. Calibration Optimality: Photon-Number Matching

Lemma 7.4 

(Scale factor calculus). Define

$$C_{2K}\left(t\right):={\displaystyle \frac{t^{-2K}}{d_{tB}}}=exp\left(-2Klogt-\frac{1}{2}\sum _{j=1}^{N}log(1-t^2{\lambda }_j^2)\right).$$

(62)

Then

$${\displaystyle \frac{d}{dt}}log{C}_{2K}\left(t\right)={\displaystyle \frac{1}{t}}\left(m_{tB}-2K\right),$$

(63)

where $m_{tB}$ is the mean photon number. In particular, ${\displaystyle \frac{d}{dt}}log{C}_{2K}\left(t\right)=0$ iff $m_{tB}=2K$, and ${\displaystyle \frac{d^2}{d{t}^2}}log{C}_{2K}\left(t\right)>0$, so $t\mapsto C_{2K}\left(t\right)$ is strictly convex on $(0,{\lambda }_1^{-1})$.

Proof. 

Differentiate:

$${\displaystyle \frac{d}{dt}}log{C}_{2K}\left(t\right)=-{\displaystyle \frac{2K}{t}}-{\displaystyle \frac{1}{2}}\sum _{j=1}^N{\displaystyle \frac{-2t{\lambda }_j^2}{1-t^2{\lambda }_j^2}}={\displaystyle \frac{1}{t}}\left(\sum _{j=1}^N{\displaystyle \frac{t^2{\lambda }_j^2}{1-t^2{\lambda }_j^2}}-2K\right)={\displaystyle \frac{m_{tB}-2K}{t}}.$$

(64)

The second derivative is the sum of a positive term ${\displaystyle \frac{2K}{t^2}}$ and strictly positive terms ${\displaystyle \frac{{\lambda }_j^2}{{(1-t^2{\lambda }_j^2)}^2}}>0$, proving strict convexity.    □

Theorem 7.5 

(Calibration Optimality). Assume ${\mathcal{R}}_{\mathrm{res}}^{\left(K\right)}\left(t\right)$ admits a representation

$${\mathcal{R}}_{\mathrm{res}}^{\left(K\right)}\left(t\right)=A{C}_{2K}\left(t\right)-B,$$

(65)

with constants $A>0$, $B\in \mathbb{R}$ independent of t. Then the unique minimizer of $t\mapsto {\mathcal{R}}_{\mathrm{res}}^{\left(K\right)}\left(t\right)$ on $(0,{\lambda }_1^{-1})$ is characterized by the photon-number matching law

$$m_{tB}=2K.$$

(66)

Proof. 

By the assumed form, ${\displaystyle \frac{d}{dt}}{\mathcal{R}}_{\mathrm{res}}^{\left(K\right)}\left(t\right)=A{\displaystyle \frac{d}{dt}}{C}_{2K}\left(t\right)$; hence critical points of ${\mathcal{R}}_{\mathrm{res}}^{\left(K\right)}$ are those of $C_{2K}$. Combine with the previous lemma and strict convexity to obtain uniqueness and $m_{tB}=2K$.    □

Remark 7.1. 

For both GBS-P and GBS-I, the leading t–dependence of the MSE/variance and thus of the residual is multiplicative by $C_{2K}\left(t\right)$, while the coefficient depends on $(a_I,B)$ but not on t. This yields the displayed form up to a harmless additive constant B (see the variance formulæ below).

7.4. Angles Quantify Conditioning: Two-Sided MSE Bounds

Let $P_{degK}$ denote the orthogonal projector onto the “blueprint” subspace spanned by the canonical basis ${\left\{e_I\right\}}_{\left|I\right|=2K}$ weighted by the coefficients $a_I$ (precisely, the range of the linear functional $\mathcal{I}$), and $P_{\mathrm{GBS}}\left(t\right)$ the projector onto the sampling feature subspace at scale t (the range of $I\mapsto Haf\left({\left(tB\right)}_I\right)$ within ${\mathcal{H}}_K$). Set

$$\gamma \left(t\right):=\parallel P_{\mathrm{GBS}}\left(t\right)P_{degK}\parallel =cos{\theta }_F\left(t\right).$$

(67)

Proposition 7.6 

(Angle–variance link). Let ${\widehat{\mu }}_K$ be either the GBS-P or GBS-I degree-K estimator from n samples at fixed t. Then there exist t–independent constants $c_{-}(a,B),c_{+}(a,B)>0$ such that

$$(1-\gamma \left(t\right))c_{-}(a,B)\le nVar\left({\widehat{\mu }}_K\right)\le (1+\gamma \left(t\right))c_{+}(a,B).$$

(68)

Equivalently, the MSE (for GBS-P, replacing Var by the leading squared bias) is sandwiched between two multiples of the split energy, with constants modulated by $(1\pm cos{\theta }_F\left(t\right))$.

Proof sketch with precise ingredients. 

Both estimators are linear statistics of the empirical law $\widehat{p}\left(t\right)\in {\mathcal{H}}_K$:

$${\widehat{\mu }}_K=\langle w,\widehat{p}\left(t\right)\rangle ,w=w^{\mathrm{P}}\mathrm{or}{w}^{\mathrm{I}}.$$

(69)

Then $nVar\left({\widehat{\mu }}_K\right)=\langle \Sigma \left(t\right)w,w\rangle$ where $\Sigma \left(t\right)$ is the covariance operator of one draw,

$$\Sigma \left(t\right)=\mathbb{E}\left[(\xi -\mathbb{E}\xi )\otimes (\xi -\mathbb{E}\xi )\right],\xi \in {\mathcal{H}}_K,\P (\xi =e_I)=p_I\left(t\right).$$

(70)

Since $\mathrm{ran}(\Sigma \left(t\right))\subset \mathrm{span}\{e_I:|I|=2K\}$ and the features live in $P_{\mathrm{GBS}}\left(t\right){\mathcal{H}}_K$, one has (up to a bounded t–independent similarity determined by $a_I$ and B)

$$\Sigma \left(t\right)\simeq P_{\mathrm{GBS}}\left(t\right)P_{degK}D{P}_{degK}{P}_{\mathrm{GBS}}\left(t\right),$$

(71)

for a positive definite diagonal D encoding the t–independent weights $I!$ and $a_I$. Two applications of standard two-subspace estimates (Friedrichs/CS-decomposition) yield

$${\lambda }_{min}\left(D\right)\parallel P_{\mathrm{GBS}}\left(t\right)P_{degK}{\parallel }^2\parallel P_{degK}w{\parallel }^2\le \langle \Sigma \left(t\right)w,w\rangle \le {\lambda }_{max}\left(D\right)\parallel P_{\mathrm{GBS}}\left(t\right)P_{degK}{\parallel }^2{\parallel P_{degK}w\parallel }^2.$$

(72)

Decompose $w=w_{\Vert }+w_{\perp }$ relative to $P_{degK}$; invoke the sharp two-subspace inequalities to bound $\parallel P_{degK}w{\parallel }^2$ above/below by $(1\pm \gamma )\parallel w_{\Vert }{\parallel }^2$ and absorb t–independent factors into $c_{\pm }(a,B)$ to conclude.    □

7.5. Lyapunov Envelope and Non-Asymptotic Rates

Theorem 7.7 

(DSFL Lyapunov for GBS). Let $(U,\mathrm{\Phi })$ denote one update step (blueprint-side map U and physical-side map Φ; counting and/or retuning). Fix $t\in (0,{\lambda }_1^{-1})$ and degree $K\in \mathbb{N}$. Assume there exists $\alpha \left(t\right)\in (0,1]$ such that, conditionally on the past,

$$\mathbb{E}\left[{\mathcal{R}}_{\mathrm{res}}^{\left(K\right)}\left(U{s}^{\left(K\right)},\mathrm{\Phi }{p}^{\left(n\right)}\left(t\right)\right)|p^{\left(n\right)}\left(t\right)\right]\le (1-\alpha \left(t\right)){\mathcal{R}}_{\mathrm{res}}^{\left(K\right)}\left(s^{\left(K\right)},p^{\left(n\right)}\left(t\right)\right).$$

(73)

Then the mean residual decays exponentially:

$$\mathbb{E}{\mathcal{R}}_{\mathrm{res}}^{\left(K\right)}\left(p^{\left(n\right)}\left(t\right)\right)\le e^{-\alpha \left(t\right)n}{\mathcal{R}}_{\mathrm{res}}^{\left(K\right)}\left(p^{\left(0\right)}\left(t\right)\right).$$

(74)

Moreover, $\alpha \left(t\right)$ is monotone improving as the principal-angle factor $\gamma \left(t\right)$ decreases, and it is maximized at the photon-number match $m_{tB}=2K$.

Proof. 

Taking total expectation in the conditional contraction yields

$$\mathbb{E}{\mathcal{R}}_{\mathrm{res}}^{\left(K\right)}\left(p^{(n+1)}\left(t\right)\right)\le (1-\alpha \left(t\right))\mathbb{E}{\mathcal{R}}_{\mathrm{res}}^{\left(K\right)}\left(p^{\left(n\right)}\left(t\right)\right).$$

(75)

Iterating and using ${(1-\alpha )}^n\le e^{-\alpha n}$ gives

$$\mathbb{E}{\mathcal{R}}_{\mathrm{res}}^{\left(K\right)}\left(p^{\left(n\right)}\left(t\right)\right)\le {(1-\alpha \left(t\right))}^n{\mathcal{R}}_{\mathrm{res}}^{\left(K\right)}\left(p^{\left(0\right)}\left(t\right)\right)\le e^{-\alpha \left(t\right)n}{\mathcal{R}}_{\mathrm{res}}^{\left(K\right)}\left(p^{\left(0\right)}\left(t\right)\right).$$

(76)

The dependence of $\alpha \left(t\right)$ on $\gamma \left(t\right)$ follows from the angle–variance control; its maximizer coincides with the minimizer of the leading residual factor, attained at $m_{tB}=2K$ by calibration optimality.    □

Corollary 7.8. 

If $n\ge \alpha {\left(t\right)}^{-1}log\left({\mathcal{R}}_{\mathrm{res}}^{\left(K\right)}\left(p^{\left(0\right)}\left(t\right)\right)/{\epsilon }^2\right)$, then $\mathbb{E}{\mathcal{R}}_{\mathrm{res}}^{\left(K\right)}\left(p^{\left(n\right)}\left(t\right)\right)\le {\epsilon }^2$. In particular, the right-hand side is minimized in t at $m_{tB}=2K$.

Remark 7.2. 

The hypothesis is a sectorial coercivity statement: in expectation, one update removes a fixed fraction of the degree-K residual. Admissibility (data–processing) and a positive principal-angle margin guarantee such a contraction. Empirically this appears as a straight line of slope $-\alpha \left(t\right)$ in semi-log MSE plots.

8. Sample-Complexity Consequences

In this section we derive, from first principles in the comparison geometry, (i) exact variance formulae for the unbiased importance estimator (GBS-I), (ii) an asymptotic mean–square error (MSE) for the probability/square–root estimator (GBS-P) under the usual sign convention, and (iii) sample–size envelopes in the $(\epsilon ,\delta )$ sense via Chebyshev/Markov bounds. Throughout we fix a degree $K\in \mathbb{N}$, a covariance B with spectrum $0<{\lambda }_N\le \cdots \le {\lambda }_1<1$, and a scale $t\in (0,{\lambda }_1^{-1})$. We also use the hafnian homogeneity

$$Haf\left({\left(tB\right)}_I\right)=t^{K}Haf\left(B_I\right),\left(\right|I|=2K),$$

(77)

and the GBS law on ${\mathcal{I}}_K:=\{I:|I|=2K\}$

$$p_I\left(t\right)={\displaystyle \frac{d_{tB}}{I!}}Haf{\left({\left(tB\right)}_I\right)}^2={\displaystyle \frac{d_{tB}{t}^{2K}}{I!}}Haf{\left(B_I\right)}^2,d_{tB}=\prod _{j=1}^N\sqrt{1-t^2{\lambda }_j^2}.$$

(78)

8.1. GBS-I (Unbiased): Exact Variance and $(\epsilon ,\delta )$ Guarantee

Define the unbiased importance estimator for ${\mu }_K^{\left(2\right)}={\sum }_{\left|I\right|=2K}{a}_{I}Haf{\left(B_I\right)}^2$ by

$${\widehat{\mu }}_K^{\mathrm{GBS}-\mathrm{I}}={\displaystyle \frac{1}{n}}\sum _{r=1}^n{Y}_r,Y\left(I\right):={\displaystyle \frac{I!a_I}{d_{tB}{t}^{2K}}},I\sim p\left(t\right).$$

(79)

Then $\mathbb{E}\left[Y\right]={\mu }_K^{\left(2\right)}$ by direct substitution. A single-line computation using $p_I\left(t\right)$ and the homogeneity gives

$$\mathbb{E}\left[Y^2\right]=\sum _{I\in {\mathcal{I}}_K}{p}_I\left(t\right)Y{\left(I\right)}^2={\displaystyle \frac{1}{d_{tB}{t}^{2K}}}\sum _{\left|I\right|=2K}{a}_I^{2}I!Haf{\left(B_I\right)}^2.$$

(80)

Proposition 8.1 

(Variance of GBS-I). For the unbiased estimator above,

$$Var\left({\widehat{\mu }}_K^{\mathrm{GBS}-\mathrm{I}}\right)={\displaystyle \frac{1}{n}}\left[{\displaystyle \frac{1}{d_{tB}{t}^{2K}}}\sum _{\left|I\right|=2K}{a}_I^{2}I!Haf{\left(B_I\right)}^2-{\left({\mu }_K^{\left(2\right)}\right)}^2\right].$$

(81)

Proof. 

Since ${\widehat{\mu }}_K^{\mathrm{GBS}-\mathrm{I}}$ is the empirical mean of i.i.d. Y with $\mathbb{E}Y={\mu }_K^{\left(2\right)}$, we have $Var\left({\widehat{\mu }}_K^{\mathrm{GBS}-\mathrm{I}}\right)=Var\left(Y\right)/n=(\mathbb{E}{Y}^2-{\left(\mathbb{E}Y\right)}^2)/n$, which after the computation of $\mathbb{E}{Y}^2$ above gives (81).    □

Corollary 8.2 

($(\epsilon ,\delta )$ sample size for GBS-I). By Chebyshev,

$$Pr\left(\left|{\widehat{\mu }}_K^{\mathrm{GBS}-\mathrm{I}}-{\mu }_K^{\left(2\right)}\right|>\epsilon \left|{\mu }_K^{\left(2\right)}\right|\right)\le {\displaystyle \frac{Var\left({\widehat{\mu }}_K^{\mathrm{GBS}-\mathrm{I}}\right)}{{\epsilon }^2{\left({\mu }_K^{\left(2\right)}\right)}^2}}.$$

(82)

Hence it suffices to take

$$n_{\mathrm{GBS}-\mathrm{I}}\ge {\displaystyle \frac{1}{{\epsilon }^2\delta }}·{\displaystyle \frac{{\displaystyle \frac{1}{d_{tB}{t}^{2K}}}{\sum }_{\left|I\right|=2K}{a}_I^{2}I!Haf{\left(B_I\right)}^2-{\left({\mu }_K^{\left(2\right)}\right)}^2}{{\left({\mu }_K^{\left(2\right)}\right)}^2}}.$$

(83)

In particular, at the photon-match $m_{tB}=2K$ (Theorem Calibration optimality), $t\mapsto d_{tB}^{-1}{t}^{-2K}$ is minimized, thus minimizing the right-hand side.

8.2. GBS-P (Bias Vanishes): Asymptotic MSE via Delta Method

For the (degree-K) probability/square–root estimator we assume the standard sign convention ${\sigma }_I=sign(Haf\left(B_I\right))$ is available.1 Define

$${\widehat{\mu }}_K^{\mathrm{GBS}-\mathrm{P}}=\sum _{\left|I\right|=2K}{a}_I{t}^{-K}\sqrt{I!}{\sigma }_I\sqrt{{\displaystyle \frac{{\widehat{p}}_I\left(t\right)}{d_{tB}}}}=g\left(\widehat{p}\left(t\right)\right),g\left(p\right)=\sum _{\left|I\right|=2K}{a}_I{t}^{-K}\sqrt{I!}{\sigma }_I\sqrt{{\displaystyle \frac{p_I}{d_{tB}}}}.$$

(84)

Since $p_I\left(t\right)=d_{tB}{t}^{2K}Haf{\left(B_I\right)}^2/I!$, one has $g\left(p\left(t\right)\right)={\mu }_K$. The multinomial delta method with covariance $Cov\left(\widehat{p}\right)={\displaystyle \frac{1}{n}}\left(\mathrm{diag}\left(p\right)-p{p}^{\top }\right)$ therefore yields

$$n\mathbb{E}\left[{\left({\widehat{\mu }}_K^{\mathrm{GBS}-\mathrm{P}}-{\mu }_K\right)}^2\right]⟶\nabla g{\left(p\right)}^{\top }\left(\mathrm{diag}\left(p\right)-p{p}^{\top }\right)\nabla g\left(p\right),(n\to \infty ).$$

(85)

A short calculation gives $\partial g/\partial p_I={\displaystyle \frac{1}{2}}{a}_I{t}^{-K}\sqrt{I!}{\sigma }_I{\left(d_{tB}^{-1}{p}_I\right)}^{-1/2}$ and hence

$$\sum _I{p}_I{\left({\displaystyle \frac{\partial g}{\partial p_I}}\right)}^2={\displaystyle \frac{1}{4d_{tB}{t}^{2K}}}\sum _{\left|I\right|=2K}{\displaystyle \frac{a_I^2}{I!}},{\left(\sum _I{p}_I{\displaystyle \frac{\partial g}{\partial p_I}}\right)}^2={\mu }_K^2.$$

(86)

We conclude:

Proposition 8.3 

(Asymptotic MSE of GBS-P). As $n\to \infty$,

$$\mathbb{E}\left[{\left|{\widehat{\mu }}_K^{\mathrm{GBS}-\mathrm{P}}-{\mu }_K\right|}^2\right]={\displaystyle \frac{1}{4n}}\left[{\displaystyle \frac{1}{d_{tB}{t}^{2K}}}\sum _{\left|I\right|=2K}{\displaystyle \frac{a_I^2}{I!}}-{\mu }_K^2\right]+o\left({\displaystyle \frac{1}{n}}\right).$$

(87)

Corollary 8.4 

($(\epsilon ,\delta )$ sample size for GBS-P). By Markov’s inequality applied to the centered square, an asymptotic sufficient sample size is

$$n_{\mathrm{GBS}-\mathrm{P}}\gtrsim {\displaystyle \frac{1}{{\epsilon }^2\delta }}·{\displaystyle \frac{{\displaystyle \frac{1}{d_{tB}{t}^{2K}}}{\sum }_{\left|I\right|=2K}{a}_I^2/I!-{\mu }_K^2}{4{\mu }_K^2}}.$$

(88)

As in Corollary 8.2, the right-hand side is minimized at the photon match $m_{tB}=2K$.

8.3. When GBS Outperforms Monte Carlo

Let ${\widehat{\mu }}_K^{\mathrm{MC}}={\displaystyle \frac{1}{n}}{\sum }_{r=1}^n{f}_K\left(X_r\right)$ be the plain Monte Carlo estimator for the degree-K contribution with $X_r\sim \mathcal{N}(0,B)$ and

$$f_K\left(x\right):=\sum _{\left|I\right|=2K}{a}_I{x}_I.$$

(89)

Then $\mathbb{E}{f}_K\left(X\right)={\mu }_K$ and a standard Wick expansion gives

$$Var\left({\widehat{\mu }}_K^{\mathrm{MC}}\right)={\displaystyle \frac{1}{n}}\left[\sum _{\begin{array}{c}\left|I\right|=2K\\ \left|J\right|=2K\end{array}}{a}_I{a}_{J}Haf\left(B_{I+J}\right)-{\mu }_K^2\right].$$

(90)

The quadratic form in (90) mixes all degree-$2K$ multi–indices pairwise via $Haf\left(B_{I+J}\right)$ and therefore grows combinatorially in K (and, under mild spectral hypotheses, exponentially in an effective dimension), whereas the GBS prefactors (81)–(87) are controlled by the single convex factor $d_{tB}^{-1}{t}^{-2K}$ after photon–number matching.

Theorem 8.5 

(Regime separation: polynomial vs exponential). Fix B with eigenvalues bounded away from 1 and coefficients $\left(a_I\right)$ on the unit sphere of ${\mathbb{R}}^{D_K}$. Then there exists $c=c\left(B\right)>0$ such that for sequences $K=K\left(N\right)$ with $K\asymp \zeta N^2$ and $\zeta >0$ one has the complexity separation

$$n_{\mathrm{GBS}}=\mathrm{poly}\left(N\right)\mathit{whereas}{n}_{\mathrm{MC}}\ge exp\left(c{N}^2\right),$$

(91)

for achieving the same $(\epsilon ,\delta )$ tolerance, provided the GBS sampler is photon–matched ($m_{tB}=2K$). For fixed small K the two complexities are of the same order, and hybrid policies (MC for small K, GBS beyond a data–driven crossover) are optimal.

Proof sketch. 

Under the matching law $m_{tB}=2K$, the GBS prefactors are minimized and remain polynomial in N (because $d_{tB}^{-1}={\prod }_j{(1-t^2{\lambda }_j^2)}^{-1/2}$ with $t<t_{max}<{\lambda }_1^{-1}$ and t bounded away from ${\lambda }_1^{-1}$ uniformly in N at fixed $\zeta$). In contrast, (90) sums $\mathrm{\Theta }\left(D_K^2\right)$ terms with $D_K\asymp N^{2K}$ and (for generic B) $Haf\left(B_{I+J}\right)$ bounded away from zero on a positive fraction of pairs, whence an exponential lower bound in $K\asymp \zeta N^2$. Careful bookkeeping of constants gives the stated separation.    □

8.4. Theory.

DPI, calibration optimality $m_{tB}=2K$, angle bounds, the exact GBS–I variance, the large-n GBS–P MSE, and the $(\epsilon ,\delta )$ envelopes are publishable as-is.

Consequences. The exact/algebraic forms (81)–(87) and the convexity identity

$${\displaystyle \frac{d}{dt}}log\left(d_{tB}^{-1}{t}^{-2K}\right)={\displaystyle \frac{m_{tB}-2K}{t}}$$

(92)

imply that photon–number matching $m_{tB}=2K$ simultaneously (i) minimizes the leading factors in the GBS variance/MSE, (ii) sharpens the principal–angle bounds (conditioning), and (iii) maximizes the Lyapunov rate in the residual contraction. These yield non–asymptotic $(\epsilon ,\delta )$ envelopes for $n_{\mathrm{GBS}}$ and the regime map against $n_{\mathrm{MC}}$ as asserted in Theorem 8.5.

9. Design Rules from DSFL (Actionable)

This section collects operational prescriptions that follow from the comparison geometry developed above and states them as verifiable statements with proofs or proof sketches. Throughout we fix a degree $K\in \mathbb{N}$, a covariance B, and a scale $t\in (0,{\lambda }_1^{-1})$. The Hilbert space, blueprint, physical law, residual and principal-angle factor are as in the previous sections.

9.1. Calibration by Photon-Number Matching

Theorem 9.1 

(Calibration rule). Let ${\mathcal{R}}_{\mathrm{res}}^{\left(K\right)}\left(t\right)$ denote the degree-K DSFL residual at scale t. Assume ${\mathcal{R}}_{\mathrm{res}}^{\left(K\right)}\left(t\right)=A{C}_{2K}\left(t\right)-B$ with $A>0$ and $C_{2K}\left(t\right)=t^{-2K}{d}_{tB}^{-1}$. Then $t^{★}$ uniquely minimizes ${\mathcal{R}}_{\mathrm{res}}^{\left(K\right)}\left(t\right)$ on $(0,{\lambda }_1^{-1})$ if and only if

$$m_{t^{★}B}=2K.$$

(93)

Moreover, at $t^{★}$ the leading factors in the variance/MSE bounds of the degree-K GBS estimators are minimized.

Proof. 

The identity ${\displaystyle \frac{d}{dt}}log{C}_{2K}\left(t\right)={\displaystyle \frac{1}{t}}(m_{tB}-2K)$ and the strict convexity of $t\mapsto C_{2K}\left(t\right)$ prove existence and uniqueness of the minimizer with stationarity condition $m_{t^{★}B}=2K$. The leading t–dependence of $nVar\left({\widehat{\mu }}_K^{\mathrm{GBS}-\mathrm{I}}\right)$ and of the large-n MSE of ${\widehat{\mu }}_K^{\mathrm{GBS}-\mathrm{P}}$ is proportional to $C_{2K}\left(t\right)$ by the variance/MSE formulæ in the sample-complexity section; hence the same minimizer applies.    □

9.2. Angle Check and Channel Selection

Define $\gamma \left(t\right)=\parallel P_{\mathrm{GBS}}\left(t\right)P_{degK}\parallel =cos{\theta }_F\left(t\right)\in [0,1)$ and let $c_{\pm }(a,B)>0$ be the t–independent constants from the angle–variance proposition.

Proposition 9.2 

(Angle gate for estimator dispersion). For the degree-K estimators at fixed t and n samples,

$$(1-\gamma \left(t\right))c_{-}(a,B)\le nVar\left({\widehat{\mu }}_K\right)\le (1+\gamma \left(t\right))c_{+}(a,B).$$

(94)

Consequently, if $\gamma \left(t\right)\le {\gamma }_{★}\in [0,1)$ then

$$nVar\left({\widehat{\mu }}_K\right)\le (1+{\gamma }_{★})c_{+}(a,B),$$

(95)

and if $\gamma \left(t\right)\ge 1-\eta$ with $\eta \in (0,1]$ then

$$nVar\left({\widehat{\mu }}_K\right)\ge \eta c_{-}(a,B).$$

(96)

Proof. 

The two-sided bound is the angle–variance link. The corollaries follow by monotonicity in $\gamma$.    □

Theorem 9.3 

(Angle threshold ⇒ GBS vs MC decision). Let $n_{\mathrm{GBS}}\left(t\right)$ and $n_{\mathrm{MC}}$ be guaranteed sample sizes (for a fixed $(\epsilon ,\delta )$) derived from Chebyshev/Markov bounds for the degree-K task. There exists a computable threshold $\overline{\gamma }=\overline{\gamma }(\epsilon ,\delta ,a,B,K)>0$ such that

$$\gamma \left(t\right)\le \overline{\gamma }⟹n_{\mathrm{GBS}}\left(t\right)\le n_{\mathrm{MC}}.$$

(97)

Conversely, if $\gamma \left(t\right)\ge 1-\eta$ with η below a computable level $\underline{\eta }(\epsilon ,\delta ,a,B,K)$, then $n_{\mathrm{GBS}}\left(t\right)\ge n_{\mathrm{MC}}$.

Proof sketch. 

Insert the angle bounds into the Chebyshev/Markov forms for $n_{\mathrm{GBS}}\left(t\right)$ (resp. the corresponding quantity for MC). The MC denominator involves the full mixed moment structure, while the GBS denominator is degree-filtered; both numerators are bounded above/below by $(1\pm \gamma )$ times t–independent constants. Solving the resulting inequality for $\gamma$ yields the thresholds $\overline{\gamma }$ and $\underline{\eta }$ (explicit in terms of $c_{\pm }(a,B)$ and the MC dispersion proxies).    □

9.2.1. Practical Surrogate

Let $X_t$ be the feature matrix whose columns are the (normalized) sampling features at scale t and $Y_K$ a basis for the degree-K blueprint range. Then $\widehat{\gamma }=\parallel X_t^{\top }{Y}_K\parallel$ satisfies $\widehat{\gamma }\ge \gamma \left(t\right)$, so $\widehat{\gamma }\le \overline{\gamma }$ is a conservative sufficient condition for choosing GBS.

9.3. Degree Decomposition and Hybrid Allocation

Let $f={\sum }_{k\ge 0}{f}_k$ be the degree decomposition and write $\mu ={\sum }_k{\mu }_k$. Suppose each degree-k component is estimated either by GBS (at its matched $t_k$) or by MC, yielding per-degree guaranteed sample sizes $n_k^{\mathrm{GBS}}$ and $n_k^{\mathrm{MC}}$ under a common $(\epsilon ,\delta )$ budget with per-degree tolerances ${\left\{{\epsilon }_k\right\}}_{k\le K}$ satisfying ${\sum }_{k\le K}{\epsilon }_k\le \epsilon$.

Proposition 9.4 

(Greedy degree-wise policy is optimal under monotone margins). Assume $n_k^{\mathrm{GBS}}$ and $n_k^{\mathrm{MC}}$ are positive and the advantage margins

$${\Delta }_k:=n_k^{\mathrm{MC}}-n_k^{\mathrm{GBS}}$$

(98)

are nondecreasing in k. Then the policy that assigns GBS to the largest k for which ${\Delta }_k>0$ and MC otherwise minimizes the total sample size ${\sum }_{k\le K}{n}_k$ among all degree-wise assignments obeying the same tolerance split.

Proof. 

Under the monotonicity hypothesis, the sequence of savings $\left\{{\Delta }_k\right\}$ is ordered; an exchange argument shows any assignment deviating from the “assign GBS to the largest positive ${\Delta }_k$” rule can be improved by swapping one GBS and one MC choice, strictly reducing total cost. Iterating yields the stated greedy optimum.    □

Corollary 9.5 

(A posteriori hybrid switch). If an a priori monotonicity check is unavailable, an empirical DSFL criterion suffices: after a warm-up with $n_0$ samples per degree,

$${\widehat{{\mathcal{R}}_{\mathrm{res}}}}_{\mathrm{GBS}}^{\left(k\right)}\le {\widehat{{\mathcal{R}}_{\mathrm{res}}}}_{\mathrm{MC}}^{\left(k\right)}⟹\mathit{use}\mathit{GBS}\mathit{at}\mathit{degree}k;\mathit{else}\mathit{use}\mathit{MC}.$$

(99)

By the Lyapunov contraction and angle bounds, this rule is asymptotically consistent.

9.4. Seam Stability Under Composition (Gluing) and Leakage Quantification

Consider a pipeline of L updates $(U_{\ell },{\mathrm{\Phi }}_{\ell })$, $\ell =1,\dots ,L$, where each $U_{\ell }$ acts on the blueprint side and each ${\mathrm{\Phi }}_{\ell }$ on the physical side (retuning, counting, linear optical blocks etc.).

Proposition 9.6 

(No inflation under gluing). If every step is DSFL-admissible, i.e. ${\mathrm{\Phi }}_{\ell }{U}_{\ell }=U_{\ell }{\mathrm{\Phi }}_{\ell }$ and $\parallel {\mathrm{\Phi }}_{\ell }\parallel \le 1$, then the composite update $(U,\mathrm{\Phi })=(U_L\cdots U_1,{\mathrm{\Phi }}_L\cdots {\mathrm{\Phi }}_1)$ satisfies

$${\mathcal{R}}_{\mathrm{res}}^{\left(K\right)}\left(U{s}^{\left(K\right)},\mathrm{\Phi }p\right)\le {\mathcal{R}}_{\mathrm{res}}^{\left(K\right)}\left(s^{\left(K\right)},p\right)\mathit{for}\mathit{all}{s}^{\left(K\right)},p.$$

(100)

Proof. 

The data–processing inequality holds at each step:

$${\mathcal{R}}_{\mathrm{res}}^{\left(K\right)}\left(U_{\ell }s,{\mathrm{\Phi }}_{\ell }p\right)\le {\mathcal{R}}_{\mathrm{res}}^{\left(K\right)}(s,p).$$

(101)

Induction on ℓ yields the composite bound.    □

Proposition 9.7 

(Angle leakage bound). Let ${\theta }_F^{(\ell )}$ be the Friedrichs angle between the degree-K blueprint subspace and the sampling feature subspace after step ℓ. If each ${\mathrm{\Phi }}_{\ell }$ is a contraction onto a feature subspace, then

$$cos{\theta }_F^{\left(L\right)}\le min\left\{1,\sum _{\ell =1}^{L}cos{\theta }_F^{(\ell )}\right\}.$$

(102)

In particular, if $cos{\theta }_F^{(\ell )}\le {\eta }_{\ell }$ with ${\sum }_{\ell }{\eta }_{\ell }<1$, then ${\theta }_F^{\left(L\right)}$ stays bounded away from 0 and the angle–variance bounds remain uniformly coercive along the pipeline.

Proof sketch. 

For orthogonal projectors $P,Q$ one has $\parallel PQ\parallel =cos{\theta }_F(P,Q)$. CS-decomposition and sub-multiplicativity give $\parallel P_{\mathrm{GBS}}^{\left(L\right)}{P}_{degK}\parallel \le \parallel P_{\mathrm{GBS}}^{\left(L\right)}{P}_{\mathrm{GBS}}^{(L-1)}\parallel \cdots \parallel P_{\mathrm{GBS}}^{\left(1\right)}{P}_{degK}\parallel .$ Each factor is bounded by the corresponding cosine; the asserted ${\ell }^1$ bound follows from the inequality $1-{\prod }_{\ell }(1-x_{\ell })\le {\sum }_{\ell }{x}_{\ell }$ applied to $x_{\ell }=1-cos{\theta }_F^{(\ell )}$.    □

9.4.1. Consequence

Under calibrated retuning at each seam (so that $m_{t_{\ell }B}=2K$), the principal angle cannot collapse in finite length pipelines formed by admissible blocks whose individual misalignments are summable; thus residuals do not inflate and the Lyapunov envelope persists.

10. Experiments

This section describes synthetic and hardware-aligned numerical studies designed to validate the DSFL predictions. We specify the experimental design, the estimators and metrics, and we state formal consistency properties that connect empirical readouts (MSE curves, angle surrogates, advantage percentages) to the theoretical quantities in the preceding sections. Throughout, degree-K quantities are taken at the photon-match $m_{tB}=2K$ unless otherwise indicated; empirical laws ${\widehat{p}}^{\left(n\right)}\left(t\right)$ are formed from n samples.

10.1. Simulated GBS and Synthetic $(B,a_I)$

10.1.1. Design.

Fix grids

$$\mathcal{N}\subset \mathbb{N},\mathcal{K}\subset \mathbb{N},\mathrm{and}\mathrm{choose}(N,K)\in \mathcal{N}\times \mathcal{K}.$$

(103)

For each $(N,K)$:

• Draw eigenvalue spectra ${\left\{{\lambda }_j\right\}}_{j=1}^N$ from prescribed classes:

  $$\left(\mathrm{flat}\right){\lambda }_j=\lambda \in (0,1),\left(\mathrm{polynomial}\right){\lambda }_j=c{j}^{-\beta },\left(\mathrm{exponential}\right){\lambda }_j=c{e}^{-\beta j},$$

  (104)

  with amplitude $c\in (0,1)$ and decay $\beta >0$ chosen so that $0<{\lambda }_N\le \cdots \le {\lambda }_1<1$.
• Draw a random Haar orthogonal $U\in O\left(N\right)$ and set $B=U\mathrm{diag}({\lambda }_1,\dots ,{\lambda }_N)U^{\top }$.
• Sample the degree-K blueprint on the unit sphere:

  $$s_K={\left(a_I\right)}_{\left|I\right|=2K}\sim \mathrm{Unif}\left({\mathbb{S}}^{D_K-1}\right),D_K=\#\{I:|I|=2K\}.$$

  (105)
• Compute the matched scale $t^{★}$ by solving $m_{tB}=2K$; equivalently, find the unique $t\in (0,{\lambda }_1^{-1})$ such that ${\sum }_{j=1}^N{t}^2{\lambda }_j^2/(1-t^2{\lambda }_j^2)=2K$.

10.1.2. Estimators and Metrics

For each configuration and for a range of sample sizes n:

• Form ${\widehat{\mu }}_K^{\mathrm{GBS}-\mathrm{I}}(n,t^{★})$ and ${\widehat{\mu }}_K^{\mathrm{GBS}-\mathrm{P}}(n,t^{★})$ from n i.i.d. synthetic GBS draws at $t^{★}$ (via the hafnian-defined pmf).
• Form the plain Monte Carlo estimator ${\widehat{\mu }}_K^{\mathrm{MC}}\left(n\right)$ using n i.i.d. $X\sim \mathcal{N}(0,B)$ and Wick reconstruction to isolate degree K if desired.
• Estimate MSEs and empirical residuals by repetitions:

  $$\widehat{\mathrm{MSE}}(•):={\displaystyle \frac{1}{R}}\sum _{r=1}^R{\left({\widehat{\mu }}^{\left(r\right)}-{\mu }_K\right)}^2,{\widehat{{\mathcal{R}}_{\mathrm{res}}}}^{\left(K\right)}\left({\widehat{p}}^{\left(n\right)}\left(t^{★}\right)\right):=\parallel {\widehat{p}}^{\left(n\right)}\left(t^{★}\right)-s_K{\parallel }_{{\mathcal{H}}_K}^2.$$

  (106)
• Report guaranteed sample sizes ${\widehat{n}}_{\mathrm{GBS}-\mathrm{I}},{\widehat{n}}_{\mathrm{GBS}-\mathrm{P}},{\widehat{n}}_{\mathrm{MC}}$ by inverting Chebyshev/Markov envelopes calibrated with empirical variances.
• Compute efficiency ratios and advantage percentages:

  $$\mathrm{Eff}.\mathrm{ratio}={\displaystyle \frac{{\widehat{n}}_{\mathrm{MC}}}{{\widehat{n}}_{\mathrm{GBS}}}},\mathrm{Adv.}\%={\displaystyle \frac{1}{T}}\sum _{t=1}^T\mathbf{1}\{{\widehat{n}}_{\mathrm{GBS}}\le {\widehat{n}}_{\mathrm{MC}}\}\times 100\%.$$

  (107)

10.1.3. Theoretical Guarantees

Proposition 10.1 

(Consistency of advantage estimators). Fix $(N,K)$ and a spectrum class with $0<{\lambda }_j<1$. Let ${\widehat{n}}_{\mathrm{GBS}}$ and ${\widehat{n}}_{\mathrm{MC}}$ be defined from empirical variances with $R\to \infty$ repetitions and $n\to \infty$ samples per repetition at matched $t^{★}$. Then

$${\widehat{n}}_{\mathrm{GBS}}\stackrel{\mathrm{a}.\mathrm{s}.}{\to }{n}_{\mathrm{GBS}}^{★},{\widehat{n}}_{\mathrm{MC}}\stackrel{\mathrm{a}.\mathrm{s}.}{\to }{n}_{\mathrm{MC}}^{★},$$

(108)

where $n^{★}$ are the population Chebyshev/Markov sample sizes. Consequently, the empirical advantage percentage converges almost surely to the true percentage under the same configuration ensemble.

Proof sketch. 

Strong laws for empirical second moments of linear statistics of i.i.d. draws imply almost-sure convergence of variance estimators. Continuity of the inversion map $(v\mapsto cv/\left({\epsilon }^2\delta \right))$ yields the stated limits; the law of large numbers over configuration draws implies convergence of frequencies.    □

10.2. Ablations Testing DSFL Predictions

10.2.1. T-Sweep (U-Shape)

For each fixed $(N,K)$ and $s_K$, sweep t in a neighborhood of $t^{★}$ and evaluate ${\widehat{{\mathcal{R}}_{\mathrm{res}}}}^{\left(K\right)}\left({\widehat{p}}^{\left(n\right)}\left(t\right)\right)$ and $\widehat{\mathrm{MSE}}\left({\widehat{\mu }}_K\left(t\right)\right)$.

Proposition 10.2 

(Empirical convexity and minimum at the match). Under the simulators above, for any fixed $(N,K)$ and sufficiently large $n,R$,

$$t\mapsto \mathbb{E}{\widehat{{\mathcal{R}}_{\mathrm{res}}}}^{\left(K\right)}\left({\widehat{p}}^{\left(n\right)}\left(t\right)\right)\mathit{and}t\mapsto \mathbb{E}\widehat{\mathrm{MSE}}\left({\widehat{\mu }}_K\left(t\right)\right)$$

(109)

are approximately convex with a unique minimum at $t=t^{★}$ (the solution of $m_{tB}=2K$). Moreover, their second derivatives with respect to $logt$ are strictly positive up to $o_n\left(1\right)$ errors.

Proof sketch. 

The residual and variance/MSE leading factors are proportional to $C_{2K}\left(t\right)=t^{-2K}{d}_{tB}^{-1}$; Lemma (scale calculus) showed ${\displaystyle \frac{d^2}{d{t}^2}}log{C}_{2K}\left(t\right)>0$ and stationarity ${\displaystyle \frac{d}{dt}}log{C}_{2K}\left(t\right)=0$ iff $m_{tB}=2K$. Empirical averages converge to their expectations by SLLN, hence the observed convex U-shape and unique minimum at $t^{★}$ up to sampling error.    □

10.2.2. Angle Surrogate vs Variance

Construct features $X_t$ and blueprint basis $Y_K$, compute $\widehat{\gamma }=\parallel X_t^{\top }{Y}_K\parallel$, and correlate $\widehat{\gamma }$ with empirical $n\widehat{Var}\left({\widehat{\mu }}_K\left(t\right)\right)$ across t.

Proposition 10.3 

(Monotone trend). For fixed $(N,K)$ and large $n,R$, the map $t\mapsto n\widehat{Var}\left({\widehat{\mu }}_K\left(t\right)\right)$ is bounded between two affine functions of $\widehat{\gamma }\left(t\right)$:

$$(1-\widehat{\gamma }\left(t\right)){\widehat{c}}_{-}(a,B)\lesssim n\widehat{Var}\left({\widehat{\mu }}_K\left(t\right)\right)\lesssim (1+\widehat{\gamma }\left(t\right)){\widehat{c}}_{+}(a,B),$$

(110)

with ${\widehat{c}}_{\pm }\to c_{\pm }$ in probability as $R\to \infty$.

Proof sketch. 

Replace the projector norm $\gamma \left(t\right)$ by the computable surrogate $\widehat{\gamma }\left(t\right)\ge \gamma \left(t\right)$; the angle–variance proposition yields the population inequalities. Empirical variances converge to their population values, and the Gram-based surrogate is continuous in sample averages of features, giving the displayed bounds up to $o_{n,R}\left(1\right)$ errors.    □

10.2.3. Hybrid Crossover

Partition the blueprint across degrees $k\le K_{max}$ and estimate the total samples needed by degree-wise GBS vs MC. Plot the crossover index $k_{\times }$ where the hybrid policy switches.

Proposition 10.4 

(Crossover index is increasing with N and decreasing with decay). Let $k_{\times }(N;\beta )$ be the smallest degree where $n_k^{\mathrm{GBS}}\le n_k^{\mathrm{MC}}$ under a polynomial (or exponential) eigen-decay parameter β. Then for sufficiently large N and fixed tolerance $(\epsilon ,\delta )$,

$${\partial }_N{k}_{\times }(N;\beta )\ge 0,{\partial }_{\beta }{k}_{\times }(N;\beta )\le 0,$$

(111)

i.e. larger systems favor GBS at lower degrees, while stronger decay (smaller effective dimension) delays the crossover.

Proof sketch. 

The MC dispersion includes mixed-moment growth with the effective dimension and increases rapidly with N when decay is weak; GBS dispersion at the matched degree depends primarily on degree-k features. Comparing the two sample-size envelopes and differentiating w.r.t. the structural parameters yields the stated monotonicities.    □

10.3. Hardware-Aligned Kernels (Optional)

10.3.1. Models

Implement time-bin/loop architectures by stochastic kernels that compose (i) squeezing maps defining $tB$, (ii) interferometer unitary draws, and (iii) photon-number resolving measurements (or threshold versions). Insert calibrated loss channels modeled as contractions on ${\mathcal{H}}_K$.

10.3.2. Stress Tests

For each pipeline, concatenate L calibrated blocks and measure ${\widehat{{\mathcal{R}}_{\mathrm{res}}}}^{\left(K\right)}$ after each seam, as well as $\widehat{\gamma }\left(t_{\ell }\right)$.

Proposition 10.5 

(No inflation under concatenation with calibrated retuning). Assume each block is DSFL-admissible and retuned to its local match $m_{t_{\ell }B}=2K$. Then the empirical residual is nonincreasing along the pipeline up to sampling error:

$$\mathbb{E}{\widehat{{\mathcal{R}}_{\mathrm{res}}}}^{\left(K\right)}\left({\widehat{p}}_{\ell }^{\left(n\right)}\right)\le \mathbb{E}{\widehat{{\mathcal{R}}_{\mathrm{res}}}}^{\left(K\right)}\left({\widehat{p}}_{\ell -1}^{\left(n\right)}\right)+o_n\left(1\right),\ell =1,\dots ,L.$$

(112)

Proof. 

Apply the DPI for each admissible block (projection/contraction and intertwining), and use convergence of empirical residuals to their expectations.    □

Proposition 10.6 

(Robustness to small losses). Let a loss channel of strength $\eta \in [0,1)$ act between calibrated blocks. Then the matched scale $t_{\ell }^{★}$ solves the perturbed equation $m_{t_{\ell }{B}^{\eta }}=2K$ and the contraction rate $\alpha \left(t_{\ell }^{★}\right)$ degrades by at most $O\left(\eta \right)$; consequently, the empirical semi-log MSE slope remains within $O\left(\eta \right)$ of the noiseless value for fixed n.

Proof sketch. 

A lossy Gaussian channel induces $B\mapsto B^{\eta }$ with $\parallel B^{\eta }-B\parallel =O\left(\eta \right)$; Davis–Kahan perturbation gives $|{\gamma }^{\eta }-\gamma |=O\left(\eta \right)$, while the scale calculus shows $t^{★}$ shifts continuously in $\eta$. The Lyapunov rate $\alpha$ is monotone in the angle margin, hence perturbs by $O\left(\eta \right)$.    □

10.3.3. Reporting

For each setting we present: (i) MSE and residual vs n (semi-log slopes), (ii) MSE vs t (U-shapes with minimum at $t^{★}$), (iii) variance vs angle surrogate curves, (iv) advantage heatmaps over $(N,K)$, and (v) pipeline residual traces across seams with/without loss.

11. Concluding Discussion

This work reframes Gaussian Boson Sampling (GBS) as a calibrated variance–reduction engine operating in a single comparison geometry. Within that geometry, the Deterministic Statistical Feedback Law (DSFL) supplies three pillars that, taken together, elevate the “folk rules’’ of GBS into theorems: (i) calibration optimality, which identifies photon–number matching $m_{tB}=2K$ as the unique minimizer of the leading variance/MSE factor $t^{-2K}{d}_{tB}^{-1}$; (ii) principal–angle conditioning, which yields sharp two–sided dispersion bounds with multiplicative constants $(1\pm cos{\theta }_F)$; and (iii) a Lyapunov envelope that converts one–step contraction into a non–asymptotic exponential error decay ${\mathrm{MSE}}_n\le C{e}^{-\alpha \left(t\right)n}$ with $\alpha \left(t\right)$ maximal at the photon match. These statements are not merely sufficient conditions; they jointly explain why GBS beats Monte Carlo (MC) in high–degree regimes, how to tune the device to realize the advantage, and what rates to expect in finite samples.

On the algorithmic side, the DSFL lens unifies unbiased (GBS–I) and biased–but–consistent (GBS–P) estimators, shows that both inherit the same convex t–dependence, and clarifies the role of degree partitioning: once the target is decomposed into degree components, GBS naturally concentrates effort where MC pays an exponential price through cross–moment couplings, while a hybrid policy remains optimal for small degrees. On the systems side, DSFL’s data–processing inequality guarantees no residual inflation under counting, retuning, or block concatenation; principal–angle control quantifies leakage across interfaces; and perturbation arguments (Davis–Kahan type) explain the graceful degradation under small hardware noise together with recovery by re–matching $m_{tB}$.

The implications are practical. First, calibration reduces to solving a strictly convex, one–dimensional equation for t, hence is robust and fast; this makes photon–number matching a pre–run, certifiable step. Second, a computable angle surrogate (e.g., a cross–Gram norm between sampling features and blueprint basis) provides an a priori gate for channel selection (GBS vs MC), enabling hybrid deployment without post–hoc tuning. Third, the Lyapunov envelope gives finite–n guarantees that can be monitored online by empirical DSFL residuals, turning semi–log MSE plots into quantitative health checks.

There are natural extensions. For non–Gaussian observables or generalized boson–sampling models (loop–hafnians, threshold detection), one must redesign the feature map and admissible updates so that the DPI persists and bias terms remain controlled; the same Hilbertian template applies. Angle surrogates can be refined—e.g. via randomized subspace embeddings—to scale to very large $(N,K)$ while preserving principal–angle fidelity. Finally, the experimental program suggested here—t–sweeps for convex U–shapes, variance–vs–angle curves, and block–wise residual traces—offers falsifiable diagnostics for hardware pipelines, and a path to controller–in–the–loop calibration where t and the interferometer are adjusted to maximize the Lyapunov rate $\alpha \left(t\right)$ subject to device constraints.

In sum, DSFL furnishes a predict–then–prove methodology for GBS: it predicts the optimal operating point ($m_{tB}=2K$), the conditioning (${\theta }_F$), and the rate ($\alpha \left(t\right)$) before data are collected; it proves that admissible updates cannot worsen the calibrated misfit; and it turns those predictions into actionable rules that persist under composition and mild noise. This synthesis clarifies when and why GBS delivers sample–complexity advantages over MC, and how to tune photonic pipelines to realize them in practice.

Appendix A.1. Standing Setup

Let $(\mathcal{H},\langle ·,·\rangle )$ be a Hilbert space with norm $\parallel x\parallel =\sqrt{\langle x,x\rangle }$, and closed subspaces $\mathsf{S},\mathsf{P}\subset \mathcal{H}$. A calibration is a pair of bounded maps

$$\mathcal{I}:\mathsf{S}\to \mathsf{P},\mathcal{J}:\mathsf{P}\to \mathsf{S},$$

(A1)

satisfying the interchangeability identities

$$\mathcal{I}\mathcal{J}={id}_{\mathsf{P}},\mathcal{J}\mathcal{I}=P_{\mathsf{S}},$$

(A2)

with $P_{\mathsf{S}}$ the orthogonal projector onto $\mathsf{S}$. For $s\in \mathsf{S}$ and $p\in \mathsf{P}$, the Residual of Sameness is

$$\mathcal{R}(s,p):={\parallel p-\mathcal{I}s\parallel }_{\mathcal{H}}^2.$$

(A3)

An admissible update is a pair $(\tilde{\mathrm{\Phi }},\mathrm{\Phi })$ with $\tilde{\mathrm{\Phi }}:\mathsf{S}\to \mathsf{S}$, $\mathrm{\Phi }:\mathsf{P}\to \mathsf{P}$ such that

$$\mathrm{\Phi }\mathcal{I}=\mathcal{I}\tilde{\mathrm{\Phi }}\mathrm{and}{\parallel \mathrm{\Phi }\parallel }_{\mathcal{H}\to \mathcal{H}}\le 1.$$

(A4)

For angle bounds, let $M,N\subset \mathcal{H}$ be closed subspaces with orthogonal projectors $P_M,P_N$, and Friedrichs cosine

$$\gamma (M,N):=\parallel P_M{P}_N\parallel \in [0,1],cos{\theta }_F(M,N)=\gamma (M,N).$$

(A5)

Appendix A.2. GBS Specialization

Let $B\in {\mathbb{R}}^{N\times N}$ be symmetric with $0<{\lambda }_N\le \cdots \le {\lambda }_1<1$, and fix a degree $K\in \mathbb{N}$. For a scale $t\in (0,{\lambda }_1^{-1})$, define

$$d_{tB}:=\prod _{j=1}^N\sqrt{1-t^2{\lambda }_j^2},C_{2K}\left(t\right):={\displaystyle \frac{t^{-2K}}{d_{tB}}},m_{tB}:=\sum _{j=1}^N{\displaystyle \frac{t^2{\lambda }_j^2}{1-t^2{\lambda }_j^2}}.$$

(A9)

Assume the degree-K residual admits the leading form ${\mathcal{R}}_{\mathrm{res}}^{\left(K\right)}\left(t\right)=A{C}_{2K}\left(t\right)-B$ with $A>0$ independent of t.

Lemma A.3 

(Scale calculus). ${\displaystyle \frac{d}{dt}log{C}_{2K}\left(t\right)=\frac{m_{tB}-2K}{t}}$, and ${\displaystyle \frac{d^2}{d{t}^2}log{C}_{2K}\left(t\right)=\frac{2K}{t^2}+\sum _{j=1}^N\frac{{\lambda }_j^2}{1-t^2{\lambda }_j^2}+\sum _{j=1}^N\frac{2t^2{\lambda }_j^4}{{(1-t^2{\lambda }_j^2)}^2}>0}$.

Proof. 

Differentiate $log{C}_{2K}\left(t\right)=-2Klogt-\frac{1}{2}{\sum }_{j}log(1-t^2{\lambda }_j^2)$ termwise. □

Theorem A.4 

(Calibration optimality). The unique minimizer $t^{★}\in (0,{\lambda }_1^{-1})$ of $t\mapsto {\mathcal{R}}_{\mathrm{res}}^{\left(K\right)}\left(t\right)$ is characterized by the photon–number match

$$m_{t^{★}B}=2K.$$

(A10)

Equivalently, $t^{★}$ minimizes the leading t–factors in the variance/MSE of the degree-K GBS estimators.

Proof. 

By Lemma A.3, critical points of $C_{2K}$ satisfy $m_{tB}=2K$ and $C_{2K}$ is strictly convex, hence the minimizer is unique. The variance/MSE leading factors for GBS-I and GBS-P are proportional to $C_{2K}\left(t\right)$; thus they share the same minimizer. □

Appendix A.3. GBS specialization (degree-K).

Let $P_{degK}$ be the blueprint projector and $P_{\mathrm{GBS}}\left(t\right)$ the sampling-feature projector at scale t, and define $\gamma \left(t\right)=\parallel P_{\mathrm{GBS}}\left(t\right)P_{degK}\parallel =cos{\theta }_F\left(t\right)$.

Theorem A.8 

(Angle–dispersion link). For the degree-K estimators from n samples at fixed t, there exist t–independent $c_{-}(a,B),c_{+}(a,B)>0$ such that

$$(1-\gamma \left(t\right))c_{-}(a,B)\le nVar\left({\widehat{\mu }}_K\right)\le (1+\gamma \left(t\right))c_{+}(a,B),$$

(A13)

with Var understood as the variance for GBS-I and as the leading asymptotic variance (or squared-bias term) for GBS-P. Equivalently, MSE is sandwiched between two multiples of the split energy, with multiplicative factors $(1\pm cos{\theta }_F\left(t\right))$.

Proof sketch. 

Write the estimator as ${\widehat{\mu }}_K=\langle w,\widehat{p}\left(t\right)\rangle$ with w fixed and $\widehat{p}\left(t\right)$ the empirical law. Then $nVar\left({\widehat{\mu }}_K\right)=\langle \Sigma \left(t\right)w,w\rangle$, where $\Sigma \left(t\right)$ is the single-sample covariance. The range of $\Sigma \left(t\right)$ lies in $P_{\mathrm{GBS}}\left(t\right)\mathcal{H}$, while w decomposes into $P_{degK}w+w_{\perp }$. Apply Proposition A.7 twice (once to control leakage into $P_{\mathrm{GBS}}\left(t\right)$, once for the split relative to $P_{degK}$) and absorb t–independent weights into $c_{\pm }(a,B)$. □

Appendix B.1. Interchangeability (Calibration)

Let $(\mathcal{H},\langle ·,·\rangle )$ be a Hilbert space carrying two closed subspaces $\mathsf{S}$ (statistical channel) and $\mathsf{P}$ (physical channel). A calibration is a pair of bounded maps

$$\mathcal{I}:\mathsf{S}\to \mathsf{P},\mathcal{J}:\mathsf{P}\to \mathsf{S}$$

(A16)

satisfying

$$\mathcal{I}\mathcal{J}={\mathrm{id}}_{\mathsf{P}},\mathcal{J}\mathcal{I}=P_{\mathsf{S}},$$

(A17)

with $P_{\mathsf{S}}$ the orthogonal projector onto $\mathsf{S}$. Thus $\mathcal{J}$ provides calibrated statistical representatives of $p\in \mathsf{P}$, $\mathcal{I}$ provides calibrated physical images of $s\in \mathsf{S}$, and composing in either channel returns the original.

Appendix B.2. One Observable: the Residual of Sameness

For $s\in \mathsf{S}$ and $p\in \mathsf{P}$,

$$\mathcal{R}(s,p):=\parallel p-\mathcal{I}\left(s\right){\parallel }_{\mathcal{H}}^2$$

(A18)

is the squared distance from p to $\mathrm{ran}\mathcal{I}$. It vanishes iff $p=\mathcal{I}\left(s\right)$ and is invariant under isometries preserving $\mathrm{ran}\mathcal{I}$.

Appendix B.3. Admissibility and the Data–Processing Inequality

A pair $(\tilde{\mathrm{\Phi }},\mathrm{\Phi })$ with $\tilde{\mathrm{\Phi }}:\mathsf{S}\to \mathsf{S}$ and $\mathrm{\Phi }:\mathsf{P}\to \mathsf{P}$ is admissible if

$$\mathrm{\Phi }\mathcal{I}=\mathcal{I}\tilde{\mathrm{\Phi }}{,\parallel \mathrm{\Phi }\parallel }_{\mathcal{H}\to \mathcal{H}}\le 1,{\parallel \tilde{\mathrm{\Phi }}\parallel }_{\mathcal{H}\to \mathcal{H}}\le 1.$$

(A19)

Then the Hilbertian data–processing inequality holds:

$$\mathcal{R}(\tilde{\mathrm{\Phi }}s,\mathrm{\Phi }p)=\parallel \mathrm{\Phi }(p-\mathcal{I}s){\parallel }_{\mathcal{H}}^2\le \parallel p-\mathcal{I}s{\parallel }_{\mathcal{H}}^2=\mathcal{R}(s,p).$$

(A20)

Hence no admissible update inflates the calibrated mismatch. In operator–algebraic models, $\mathcal{J}$ is an orthogonal conditional expectation, and (34) follows from Kadison–Schwarz.

Appendix B.4. Propagation with a Gap: the DSFL Lyapunov Envelope

Assume a residual identity

$${\displaystyle \frac{d}{dt}}\mathcal{R}\left(t\right)=-2{\langle Ke,e\rangle }_{\mathcal{H}}+2{\langle e,g\rangle }_{\mathcal{H}},e:=p-\mathcal{I}s,K=K^{*}⪰0,$$

(A21)

and bounds

$${\langle Ke,e\rangle }_{\mathcal{H}}\ge {\kappa \parallel e\parallel }_{\mathcal{H}}^2{,|\langle e,g\rangle }_{\mathcal{H}}{|\le \epsilon \parallel e\parallel }_{\mathcal{H}}^2,\kappa >\epsilon \ge 0.$$

(A22)

Then

$${\displaystyle \frac{d}{dt}}\mathcal{R}\left(t\right)\le -\alpha \mathcal{R}\left(t\right),\alpha :=2(\kappa -\epsilon )>0,\Rightarrow \mathcal{R}\left(t\right)\le e^{-\alpha (t-t_0)}\mathcal{R}\left(t_0\right).$$

(A23)

Thus $log\mathcal{R}$ decays at a single rate $\alpha$ fixed by the sector’s coercivity margin.

Appendix B.5. Operational Arrow (Accumulated Improvement)

Define

$$S_R\left(t\right):=-log{\displaystyle \frac{\mathcal{R}\left(t\right)}{\mathcal{R}\left(0\right)}}={\int }_0^t\alpha \left(\tau \right)d\tau .$$

(A24)

It is nondecreasing and quantifies total calibrated mismatch removed, independently of the choice of clock.

Appendix B.6. Minimal Checklist for Applications

• Fix $\mathcal{H}$, $\mathsf{S}$, $\mathsf{P}$; choose $(\mathcal{I},\mathcal{J})$ with (30).
• Verify (A19) for each step; (34) then guarantees monotonicity.
• If dynamics are available, identify $(\kappa ,\epsilon )$ in (A22) to obtain $\alpha$ in (A23).

Appendix B.7. Generic Two–Channel Application Template

Let $s\left(t\right)\in \mathsf{S}$ and $p\left(t\right)\in \mathsf{P}$ with calibration $\mathcal{I}:\mathsf{S}\to \mathsf{P}$. Set $e\left(t\right):=p\left(t\right)-\mathcal{I}s\left(t\right)$, $\mathcal{R}\left(t\right):={\parallel e\left(t\right)\parallel }_{\mathcal{H}}^2$. If

$$\dot{\mathcal{R}}\left(t\right)=-2{\langle Ke,e\rangle }_{\mathcal{H}}+2{\langle e,g\rangle }_{\mathcal{H}},{\langle Ke,e\rangle }_{\mathcal{H}}\ge {\kappa \parallel e\parallel }_{\mathcal{H}}^2{,|\langle e,g\rangle }_{\mathcal{H}}{|\le \epsilon \parallel e\parallel }_{\mathcal{H}}^2,$$

(A25)

then $\dot{\mathcal{R}}\le -2(\kappa -\epsilon )\mathcal{R}$ and $\mathcal{R}\left(t\right)\le e^{-2(\kappa -\epsilon )(t-t_0)}\mathcal{R}\left(t_0\right)$. This dynamic envelope pairs with the static monotonicity $\mathcal{R}(\tilde{\mathrm{\Phi }}s,\mathrm{\Phi }p)\le \mathcal{R}(s,p)$ whenever $(\tilde{\mathrm{\Phi }},\mathrm{\Phi })$ is admissible.

Appendix B.8. Terminology and Mechanism: sDoF/pDoF, DSFL, and the Dual Feedback Loop

Appendix B.8.1. What are sDoF and pDoF?

The statistical degrees of freedom (sDoF) are the blueprint variables $s\in \mathsf{S}$ that encode the model, prior, or control law; the physical degrees of freedom (pDoF) are the measured/realized responses $p\in \mathsf{P}$. Both live in the same comparison Hilbert space $(\mathcal{H},\langle ·,·\rangle )$ via a calibration $(\mathcal{I},\mathcal{J})$ obeying

$$\mathcal{I}\mathcal{J}={\mathrm{id}}_{\mathsf{P}},\mathcal{J}\mathcal{I}=P_{\mathsf{S}}.$$

(A26)

Equation (A26) ensures two-way coherence: $\mathcal{J}$ reads a calibrated blueprint from a physical state, and $\mathcal{I}$ writes a calibrated physical image of a statistical state. The Residual of Sameness

$$\mathcal{R}(s,p)={\parallel p-\mathcal{I}s\parallel }_{\mathcal{H}}^2$$

(A27)

measures the mismatch between pDoF and the calibrated sDoF within $\mathcal{H}$.

Appendix B.8.2. Why “DSFL”

Deterministic Statistical Feedback Law (DSFL) is shortened to DSFL because updates are deterministic (linear/contractive), the driver is statistical (sDoF), the mechanism is feedback (via the residual $e:=p-\mathcal{I}s$), and it functions as a law (one inequality governs all steps). For any admissible update $(\tilde{\mathrm{\Phi }},\mathrm{\Phi })$ ($\mathrm{\Phi }\mathcal{I}=\mathcal{I}\tilde{\mathrm{\Phi }}$, $\parallel \mathrm{\Phi }\parallel ,\parallel \tilde{\mathrm{\Phi }}\parallel \le 1$),

$$\mathcal{R}(\tilde{\mathrm{\Phi }}s,\mathrm{\Phi }p)=\parallel \mathrm{\Phi }(p-\mathcal{I}s){\parallel }_{\mathcal{H}}^2\le {\parallel p-\mathcal{I}s\parallel }_{\mathcal{H}}^2=\mathcal{R}(s,p),$$

(A28)

so the residual is $L^2$–monotone.

Appendix A.8.3. Dual Feedback Loop (Discrete Time)

There are two tightly coupled feedback paths in the same geometry:

$$\text{S-loop}:s_{k+1}=\tilde{\mathrm{\Phi }}{s}_k,\text{P-loop}:p_{k+1}=\mathrm{\Phi }{p}_k,e_{k+1}=p_{k+1}-\mathcal{I}{s}_{k+1},$$

(A29)

with $\mathrm{\Phi }\mathcal{I}=\mathcal{I}\tilde{\mathrm{\Phi }}$ and $\parallel \mathrm{\Phi }\parallel ,\parallel \tilde{\mathrm{\Phi }}\parallel \le 1$. By (A28), $\mathcal{R}(s_{k+1},p_{k+1})\le \mathcal{R}(s_k,p_k)$ sample–wise. If, in addition, a Lyapunov identity with a gap holds, $\dot{\mathcal{R}}\le -\alpha \mathcal{R}$, then ${\mathcal{R}}_k\le C{e}^{-\alpha k}$ (exponential envelope).

Appendix B.8.4. Flow form (Continuous Time)

Writing $e\left(t\right)=p\left(t\right)-\mathcal{I}s\left(t\right)$, assume

$${\displaystyle \frac{d}{dt}}\mathcal{R}\left(t\right)=-2{\langle Ke,e\rangle }_{\mathcal{H}}+2{\langle e,g\rangle }_{\mathcal{H}},\langle Ke,e\rangle \ge {\kappa \parallel e\parallel }^2,\left|\langle e,g\rangle \right|\le {\epsilon \parallel e\parallel }^2,$$

(A30)

with $\kappa >\epsilon \ge 0$. Then ${\displaystyle \frac{d}{dt}}\mathcal{R}\le -2(\kappa -\epsilon )\mathcal{R}$ and $\mathcal{R}\left(t\right)\le e^{-2(\kappa -\epsilon )(t-t_0)}\mathcal{R}\left(t_0\right)$. The accumulated drop $S_R\left(t\right)=-log{\displaystyle \frac{\mathcal{R}\left(t\right)}{\mathcal{R}\left(0\right)}}={\int }_0^{t}2(\kappa -\epsilon )d\tau$ is nondecreasing and quantifies the total calibrated mismatch removed.

Appendix B.8.5. Practical Interpretation

The P-loop evolves pDoF through an $L^2$–contractive map that respects calibration, reducing e on the physical side; the S-loop refines sDoF so that $\mathcal{I}s$ tracks the evolving p. Because both loops live in $\mathcal{H}$ and intertwine $(\mathcal{I},\mathcal{J})$, they cannot fabricate residual—only dissipate or transport it. In applications with slow drift, one may update the calibration $(\mathcal{I},\mathcal{J})\mapsto ({\mathcal{I}}^{\left(m\right)},\mathcal{J}\left(m\right))$ on a coarser timescale as long as ${\mathcal{I}}^{\left(m\right)}\mathcal{J}\left(m\right)={\mathrm{id}}_{\mathsf{P}}$, $\mathcal{J}\left(m\right){\mathcal{I}}^{\left(m\right)}=P_{\mathsf{S}}$; between updates the inner dual loop remains admissible and residual–monotone.

Appendix B.8.6. Data–Processing Inequality (DPI) in the DSFL Geometry

Claim. Admissible updates (intertwining $\mathrm{\Phi }\mathcal{I}=\mathcal{I}\tilde{\mathrm{\Phi }}$ and nonexpansive ${\parallel \mathrm{\Phi }\parallel }_{\mathcal{H}\to \mathcal{H}}\le 1$) cannot inflate the residual.

Status. Proved in the text (Theorem “DPI”) by a one–line Hilbert–space contraction; stability under composition is immediate.

Appendix B.8.7. Calibration Optimality m_tB = 2K

Claim. The unique minimizer of the leading variance/MSE factor is the photon–number match.

Status. Proved by strict convexity on $(0,{\lambda }_1^{-1})$ of

$$log{C}_{2K}\left(t\right)=-2Klogt-\frac{1}{2}\sum _{j}log\left(1-t^2{\lambda }_j^2\right),{\displaystyle \frac{d}{dt}}log{C}_{2K}\left(t\right)={\displaystyle \frac{m_{tB}-2K}{t}}.$$

(A31)

Scope. Optimal for any quantity whose t–dependence factors through $C_{2K}\left(t\right)$: the leading terms of the GBS-I variance and GBS-P MSE, and the DSFL residual representation ${\mathcal{R}}_{\mathrm{res}}^{\left(K\right)}\left(t\right)=A{C}_{2K}\left(t\right)-B$ (when that representation holds).

Appendix B.8.8. Angle (Friedrichs) Bounds

Claim. Two–sided MSE/variance control with factors $(1\pm cos{\theta }_F)$.

Status. Proved via the two–subspace inequalities and $\parallel P_M{P}_N\parallel =cos{\theta }_F$; explicit derivations included.

Appendix B.8.9. Exact Variance for GBS-I; Asymptotic MSE for GBS-P

GBS-I. Exact variance formula proved from the pmf and hafnian homogeneity.

GBS-P. Large-n leading term proved by the multinomial delta method (gradient and covariance computed).

Scope/assumptions. For GBS-P we assume the usual sign-known convention and interior probabilities/smoothness so the delta method applies.

Appendix B.8.10. Chebyshev/Markov (ε,δ) Guarantees

Claim. Invert variance/MSE to obtain guaranteed sample sizes.

Status. Proved once the variance/MSE statements hold (and $\mu \ne 0$), by direct Chebyshev/Markov bounds.

Appendix B.8.11. Lyapunov Envelopes (Discrete/Continuous)

Claim. If a one–step conditional contraction holds, the expected residual decays exponentially.

Status. Proved as an implication via discrete (and continuous) Grönwall.

Scope. The contraction is the hypothesis; the envelope follows rigorously.

Appendix B.8.12. Residual Representation

Where used. To map calibration optimality to residual minimization.

Status. Provable for degree-K laws and the two estimators, since the only nontrivial t–dependence enters through $t^{-2K}$ and $d_{tB}^{-1}$. One factors t out of (i) the multinomial covariance and (ii) the linear readout; a short lemma closes the argument.

Appendix B.8.13. “Angle ⇒ rate” Linkage

Claim. The Lyapunov rate $\alpha \left(t\right)$ improves as $\gamma \left(t\right)=cos{\theta }_F\left(t\right)$ decreases.

Status. Provable once the conditional contraction constant is written as a monotone functional of the split energy controlled by $(1\pm \gamma )$. With the projector bounds in hand, a brief comparison lemma formalizes this.

Appendix B.8.14. Composition/Gluing Stability and Small-Noise Robustness

Claim. No inflation under concatenation; small perturbations degrade angles/rates by $O(\parallel \Delta \parallel )$.

Status. Concatenation: already follows from DPI under composition (proved). Robustness: Davis–Kahan/Wedin gives

$$|sin{\theta }_F(B+\Delta )-sin{\theta }_F\left(B\right)|=O\left({\displaystyle \frac{\parallel \Delta \parallel }{\mathrm{gap}}}\right),$$

(A32)

hence $\gamma$ perturbs by $O(\parallel \Delta \parallel )$ and rates follow from the previous item.

Appendix B.8.15. Asymptotic Normality (CLT) for GBS-P

Claim. Beyond the leading MSE constant, a full CLT holds.

Status. Provable by standard smooth-functionals-of-the-multinomial arguments under mild regularity (already used for the delta-method step).

Appendix B.8.16. Open, Nonempty Regions of Advantage; Percentage-of-Advantage

Content (Andersen–Shan 2025). Large-n leading term for GBS-P (matching ours), improved tuning via $m_{tB}=2K$, $(\epsilon ,\delta )$ sample-size envelopes for both estimators, existence of nonempty open subsets with advantage over MC, and percentage-of-advantage values over $(a,B)$ under Haar–Vandermonde measures.

Status. We rely on these theorems and cite them explicitly.

Appendix B.8.17. Principal-angle surrogate $\hat{\gamma }={\parallel X_t^{\top }{Y}_K\parallel }_2$ as a certified upper bound.

Claim. The computable Gram surrogate upper-bounds the true Friedrichs cosine.

Status. Heuristic but tighten-able: if $X_t$ spans the sampling feature subspace and $Y_K$ spans the blueprint subspace (properly normalized), then $\parallel X_t^{\top }{Y}_K\parallel \uparrow \parallel P_{\mathrm{GBS}}\left(t\right)P_{degK}\parallel$. A short approximation lemma with sampling-error bounds can make this rigorous.

Appendix B.8.18. Regime Separation (Poly vs. Exp).

Claim. For $K\sim \zeta N^2$, $n_{\mathrm{GBS}}=\mathrm{poly}\left(N\right)$ while $n_{\mathrm{MC}}=exp\left(\mathrm{\Theta }\left(N\right)\right)$.

Status. Proof sketch given; a fully rigorous version needs fixed ensembles (eigenvalues in $[c,1-\epsilon ]$; coefficients on the unit sphere) and a positive-mass lower bound for the MC mixed-hafnian quadratic form. With these in place, the exponential lower bound follows by counting/Stirling estimates (in line with Andersen–Shan open-set results).

Appendix B.8.19. One-Step Conditional Contraction (Hypothesis in Lyapunov Theorem)

Claim. Each update removes a fixed fraction of the residual in expectation.

Status. Currently an assumption (empirically supported by semi-log linearity). A device-level proof would require an explicit Markov kernel and a uniform spectral gap; feasible but model-specific.

Appendix B.8.20. One-Step Conditional Contraction (Hypothesis in Lyapunov Theorem)

Claim. Each update removes a fixed fraction of the residual in expectation.

Status. Currently an assumption (empirically supported by semi-log linearity). A device-level proof would require an explicit Markov kernel and a uniform spectral gap; feasible but model-specific.

Lemma B.2 

(Rate improves as angle opens). Let ${\mathcal{R}}_{\mathrm{res}}^{\left(K\right)}\left(t\right)={\parallel p-\mathcal{I}s\parallel }^2$ be the degree-K residual at scale t and let $(U,\mathrm{\Phi })$ be an admissible one-step update (intertwining and $\parallel \mathrm{\Phi }\parallel \le 1$). Assume the conditional one-step contraction hypothesis of Theorem A.9: there exists $\alpha \left(t\right)\in (0,1]$ such that

$$\mathbb{E}\left[{\mathcal{R}}_{\mathrm{res}}^{\left(K\right)}\left(Us,\mathrm{\Phi }p\right)|s,p\right]\le \left(1-\alpha \left(t\right)\right){\mathcal{R}}_{\mathrm{res}}^{\left(K\right)}\left(t\right).$$

(A41)

Then there exist t–independent constants $0<c_{-}\le c_{+}<\infty$ (depending only on the fixed weights $a_I$ and B through bounded positive operators) such that

$$c_{-}(1-\gamma \left(t\right))\le \alpha \left(t\right)\le c_{+}(1-\gamma \left(t\right)),\gamma \left(t\right)=∥P_{\mathrm{GBS}}\left(t\right)P_{degK}∥=cos{\theta }_F\left(t\right).$$

(A42)

In particular, $\alpha \left(t\right)$ increases as $\gamma \left(t\right)$ decreases and is maximized at the photon–number match $m_{tB}=2K$.

Proof. 

Write the residual as ${\mathcal{R}}_{\mathrm{res}}^{\left(K\right)}\left(t\right)={\parallel e\parallel }^2$ with $e:=p-\mathcal{I}s\in {\mathcal{H}}_K$. One step gives

$$e^{+}:=\mathrm{\Phi }p-\mathcal{I}Us=\mathrm{\Phi }(p-\mathcal{I}s)+\underset{=0\mathrm{by}\mathrm{intertwining}}{\underbrace{\mathrm{\Phi }\mathcal{I}s-\mathcal{I}Us}}=\mathrm{\Phi }e,$$

(A43)

hence ${{\mathcal{R}}_{\mathrm{res}}^{\left(K\right)}}^{+}=\parallel e^{+}{\parallel }^2={\parallel \mathrm{\Phi }e\parallel }^2$. Taking conditional expectation and comparing to ${\parallel e\parallel }^2$,

$$\alpha \left(t\right)=1-{\displaystyle \frac{{\mathbb{E}\parallel \mathrm{\Phi }e\parallel }^2}{{\parallel e\parallel }^2}}.$$

(A44)

Decompose $e=e_{\Vert }+e_{\perp }$ with respect to the blueprint/sampling pair,

$$e_{\Vert }:=P_{degK}e,e_{\perp }:=(I-P_{degK})e,$$

(A45)

and note that $\mathrm{\Phi }$ acts through the sampling feature subspace (by admissibility) so that the “through–subspace’’ contraction applies:

$$\parallel \mathrm{\Phi }e\parallel =\parallel P_{\mathrm{GBS}}\left(t\right)\mathrm{\Phi }{P}_{degK}{e}_{\Vert }+P_{\mathrm{GBS}}\left(t\right)\mathrm{\Phi }{e}_{\perp }\parallel \le \parallel \mathrm{\Phi }\parallel \parallel P_{\mathrm{GBS}}\left(t\right)P_{degK}{e}_{\Vert }\parallel +\parallel \mathrm{\Phi }\parallel \parallel P_{\mathrm{GBS}}\left(t\right)e_{\perp }\parallel .$$

(A46)

Using $\parallel \mathrm{\Phi }\parallel \le 1$ and the two–subspace bounds (Proposition A.7),

$$\parallel P_{\mathrm{GBS}}\left(t\right)P_{degK}{e}_{\Vert }\parallel \le \gamma \left(t\right)\parallel e_{\Vert }\parallel ,\parallel P_{\mathrm{GBS}}\left(t\right)e_{\perp }\parallel \le \parallel e_{\perp }\parallel .$$

(A47)

Squaring and taking expectation (conditionally on e) yields

$${\mathbb{E}\parallel \mathrm{\Phi }e\parallel }^2\le \gamma {\left(t\right)}^2\parallel e_{\Vert }{\parallel }^2+\parallel e_{\perp }{\parallel }^2\le \left(1-(1-\gamma {\left(t\right)}^2)\frac{\parallel e_{\Vert }{\parallel }^2}{{\parallel e\parallel }^2}\right){\parallel e\parallel }^2.$$

(A48)

By the two–subspace lower bound,

$$(1-\gamma \left(t\right))\left(\parallel e_{\Vert }{\parallel }^2+{\parallel e_{\perp }\parallel }^2\right)\le \parallel e_{\Vert }+e_{\perp }{\parallel }^2={\parallel e\parallel }^2\Rightarrow {\displaystyle \frac{\parallel e_{\Vert }{\parallel }^2}{{\parallel e\parallel }^2}}\ge {\displaystyle \frac{1-\gamma \left(t\right)}{2}},$$

(A49)

where the constant $\frac{1}{2}$ is an absolute (non-optimized) choice; any fixed $c\in (0,1)$ would do after renormalizing constants. Hence

$${\displaystyle \frac{{\mathbb{E}\parallel \mathrm{\Phi }e\parallel }^2}{{\parallel e\parallel }^2}}\le 1-\underset{\asymp (1-\gamma (t\left)\right)}{\underbrace{\frac{1}{2}(1-\gamma {\left(t\right)}^2)(1-\gamma \left(t\right))}},$$

(A50)

which implies an upper bound

$$\alpha \left(t\right)\ge c_{-}(1-\gamma \left(t\right)),$$

(A51)

for some $c_{-}>0$ independent of t (absorbing the numerical constants and the fixed weights). A matching upper bound follows by reversing the roles of the projectors in the two–subspace inequalities:

$$\parallel P_{\mathrm{GBS}}\left(t\right)P_{degK}{e}_{\Vert }\parallel \ge \sqrt{1-\gamma {\left(t\right)}^2}\parallel e_{\Vert }{\parallel \Rightarrow \mathbb{E}\parallel \mathrm{\Phi }e\parallel }^2\ge \left(1-(1-\gamma {\left(t\right)}^2)\frac{\parallel e_{\Vert }{\parallel }^2}{{\parallel e\parallel }^2}\right){\parallel e\parallel }^2,$$

(A52)

and using $\frac{\parallel e_{\Vert }{\parallel }^2}{{\parallel e\parallel }^2}\le \frac{1}{1-\gamma \left(t\right)}$ (another consequence of Proposition A.7) gives

$$\alpha \left(t\right)\le c_{+}(1-\gamma \left(t\right))$$

(A53)

for some $c_{+}<\infty$ independent of t. The monotonicity of $\alpha \left(t\right)$ in $\gamma \left(t\right)$ is immediate, and maximization at $m_{tB}=2K$ follows since $\gamma \left(t\right)$ is minimized (best alignment) at the photon–match by the calibration optimality and angle–conditioning discussion. □

Appendix B.9. Bottom Line

Rock solid now: DPI; calibration optimality $m_{tB}=2K$; angle bounds; exact GBS-I variance; large-n GBS-P MSE; $(\epsilon ,\delta )$ inversions; Lyapunov envelope (conditional on a one-step contraction); composition stability; small-noise robustness (Davis–Kahan).

Provable with short add-ons: residual factorization $A{C}_{2K}-B$; monotone “angle ⇒ rate” linkage; CLT for GBS-P; certified Gram surrogate bound.

Backed by the literature: nonempty open regions with advantage and percentage-of-advantage computations (Andersen–Shan 2025).

Still to complete: a model-independent proof of the one-step contraction and a fully quantified “poly vs. exp” lower bound with explicit constants under clear ensembles.