Can pure bosonic Gaussian states be agnostically tomographed efficiently?

Agnostic tomography targets the best Gaussian fit, not the state

The paper studies agnostic tomography of pure bosonic Gaussian states. Given copies of an arbitrary n-mode bosonic state ρ, the algorithm must output a pure Gaussian state |ψ̂⟩ whose infidelity with ρ is at most opt + ε. Here opt = 1 − sup_{|ψ⟩∈G} ⟨ψ|ρ|ψ⟩ is the minimum infidelity achievable by any pure Gaussian state, and G is the class of pure bosonic Gaussian states. Infidelity is defined as d_IF(ρ, σ) ≡ 1 − F(ρ, σ), with fidelity F(ρ, σ) ≡ (tr √(√ρ σ √ρ))².

This is an agnostic problem because the algorithm is not promised that ρ is Gaussian. It must compete against the best Gaussian approximation that exists for that particular ρ, whatever ρ happens to be. The guarantee is stated as d_IF(ρ, |ψ̂⟩) ≤ opt + ε, achieved with probability at least 0.99.

The distinction matters for how the result should be read. A small opt means ρ is close to Gaussian and the output is a good description of the state. A large opt means ρ is far from Gaussian, and the output is the closest Gaussian skeleton rather than a faithful reconstruction. The paper treats both regimes separately because they require different algorithms and yield different complexity.

High fidelity is strongly polynomial, low fidelity is not

In the high-fidelity regime, the paper assumes opt is below some universal constant. Under that assumption, there is an algorithm that takes poly(n, ε⁻¹, log log E) copies and uses only Gaussian unitaries, heterodyne measurements, and photon counting, and outputs |ψ̂⟩ ∈ G with d_IF(ρ, |ψ̂⟩) ≤ opt + ε, succeeding with probability at least 0.99. Here E is the energy of some closest state in G. The runtime and copy complexity are strongly polynomial in n, 1/ε and log log E. The authors describe this as the first robust Gaussian state learning guarantee for any value of opt bounded away from 0.

For arbitrary opt, including the low-fidelity regime where the best Gaussian fidelity is an arbitrary constant bounded away from one, the protocol uses (n+1) poly(1/ε) poly(1 + log log E) copies and runtime. The guarantee d_IF(ρ, |ψ̂⟩) ≤ opt + ε still holds with probability at least 0.99, but the dependence on the target precision ε is exponential rather than polynomial.

The exponential dependence is not an artifact of the analysis. The paper proves that agnostic tomography of pure Gaussian states is at least as hard as finding maximal cliques in graphs, so unless NP ⊆ BQP there is no algorithm achieving runtime poly(n, 1/ε). The two regimes therefore differ in a way readers should keep separate: the high-fidelity regime is efficient in both n and 1/ε, while the general regime is efficient in n but necessarily exponential in 1/ε.

The first tolerant Gaussianity tester for bosonic states

A direct corollary of the low-fidelity protocol is a tolerant tester for Gaussianity. The tester distinguishes states for which opt > c + ε from states for which opt < c − ε, for any threshold c ∈ (0, 1). This is the first guarantee of its kind for bosonic states.

The comparison the paper draws is with recent Gaussianity testing algorithms, which only tolerate small constant-infidelity deviations from Gaussianity. Those protocols can detect when a bosonic state is close to or far from being Gaussian, but they cannot distinguish fidelity above or below a specified threshold, even a threshold as permissive as 1/2. The new tester works for any threshold c in the open unit interval, which is what makes it tolerant rather than a fixed-threshold test.

The tester inherits the complexity of the underlying protocol, so its ε-dependence is exponential in the general regime. The relevant reader-facing point is the threshold structure: the guarantee is about separating opt above c + ε from opt below c − ε, not about estimating opt to additive error ε.

Gaussian measurements give a warm start, non-Gaussian ones refine it

Both protocols follow a shared two-step paradigm. First, the algorithm iteratively uses general Gaussian measurements combined with techniques from classical robust statistics to obtain a good warm-start estimate. Heterodyne measurements are the key primitive: when applied to a Gaussian state, the measurement output is a sample from a classical high-dimensional Gaussian whose parameters correspond directly to the parameters of the underlying state and the specific Gaussian measurement applied. Applied to a possibly non-Gaussian ρ, the outcome is a draw from a classical distribution close in total variation distance to a Gaussian whose parameters would recover the best bosonic Gaussian approximation to ρ.

This reduction places the problem inside the oblivious contamination model of robust statistics, and the agnostic problem inherits the difficulties of that setting. In the high-fidelity regime, the reduction becomes the standard problem of robustly learning a Gaussian, for which polynomial-time algorithms with strong error guarantees are known. Iteratively combining robust Gaussian estimators with techniques and bounds from prior work yields a pure Gaussian state whose infidelity with ρ is Õ(opt), using only Gaussian measurements and poly(n, log log E) copies.

The second step leverages non-Gaussian measurements, specifically photon counting, to refine the warm start using convex and non-convex optimization methods. The warm start is used to rotate the state so that, in the new basis, the best Gaussian approximation approximately lives in a low photon number subspace. Finding the closest Gaussian state then becomes a tensor optimization problem in that subspace.

List decoding handles the low-fidelity regime

The low-fidelity regime is harder because the fidelity α = 1 − opt of ρ with its best Gaussian approximation π* may be small. In that setting it is information-theoretically impossible to output a single good estimate. The paper instead uses list-decodable learning, applying a list-decodable Gaussian learner to the heterodyne outcomes to obtain a list, of size depending only on α, of candidate means and covariances. At least one entry approximates the parameters of π* up to factors depending only on α.

Such a coarse estimate is not by itself a good Gaussian approximation, but it is good enough to unsqueeze the state: it identifies a rotated basis in which the true state has smaller energy. The authors develop a list-decodable variant of iterative unsqueezing that follows every candidate in the list through the corresponding change of basis, then prunes the resulting bases using a photon-number-based notion of distance between bases, so that only O(1/α) bases survive each round.

After O(log log E) rounds, this produces a short list of rotated bases, and in at least one of them π* has mean photon number at most poly(1/α), independent of both n and E. The pruning step is what keeps the list from growing across rounds; the paper bounds the number of surviving bases through a packing argument and shows that useful bases survive a mass test when the cutoff photon number k satisfies k ≥ C₀ h/α with C₀ = 64.

Non-Gaussian measurements are provably necessary

The protocols rely on photon counting, a non-Gaussian measurement. The paper shows this reliance is not incidental. In the high-fidelity regime, there is a time-efficient protocol that uses only Gaussian measurements and achieves infidelity O(opt) + ε using Õ(n²/ε) copies. However, there exist settings of opt such that any protocol in this regime using only Gaussian measurements cannot achieve infidelity opt + ε without using exp(Ω(n^{1/3})) copies.

This is a separation between Gaussian and non-Gaussian measurements for this task. The Gaussian-measurement protocol reaches O(opt) + ε, which is a semi-agnostic guarantee, while the full protocol reaches opt + ε, the strong agnostic guarantee. The exponential lower bound shows the gap cannot be closed by a better Gaussian-measurement algorithm.

The authors also note that the strong agnostic guarantees they obtain are provably superior to what is possible for robustly estimating classical Gaussians. The comparison is between the bosonic problem and its classical analogue, and it runs in the direction of the quantum problem admitting stronger guarantees than the classical one.

Where the guarantee stops and what to weigh

The guarantee is stated for pure bosonic Gaussian states as the output class. The algorithm outputs |ψ̂⟩ ∈ G, a pure Gaussian state, and the error is measured as infidelity against ρ. The bound d_IF(ρ, |ψ̂⟩) ≤ opt + ε is relative to opt, the best infidelity achievable by any pure Gaussian state for that same ρ, so the guarantee is about matching the best possible Gaussian fit rather than about absolute reconstruction accuracy.

The high-fidelity result requires opt to be below a universal constant; the strongly polynomial dependence on 1/ε is not claimed outside that regime. The general result holds for arbitrary opt but pays an exponential dependence on 1/ε, and the hardness reduction shows this is unavoidable unless NP ⊆ BQP. The tolerant tester inherits the same ε-dependence.

The measurement model is also part of the scope. The protocols use Gaussian unitaries, heterodyne measurements, and photon counting. The separation result shows that dropping photon counting and using only Gaussian measurements costs exp(Ω(n^{1/3})) copies in the settings constructed in the paper. Readers comparing this work with Gaussian-measurement-only approaches should treat that copy-complexity gap, not just the infidelity target, as the operative difference.

The priority is the opt + ε target, not just polynomial n

The most decision-relevant feature of this paper is the target it achieves: infidelity opt + ε, where opt is the best any pure Gaussian state can do for the given ρ. Earlier Gaussianity testing algorithms tolerate only small constant-infidelity deviations, and Gaussian-measurement-only protocols reach O(opt) + ε rather than opt + ε. The step from O(opt) + ε to opt + ε is what the non-Gaussian measurements buy, and the exp(Ω(n^{1/3})) lower bound shows the step is not free.

The second feature is the split between the two regimes. If opt is below the universal constant, the reader can expect runtime and copy complexity strongly polynomial in n, 1/ε and log log E. If opt is arbitrary, the reader should expect (n+1) poly(1/ε) poly(1 + log log E) copies and runtime, with the exponential dependence on 1/ε being provably necessary unless NP ⊆ BQP.

For a reader deciding whether this approach applies to a given bosonic state, the operative question is where that state's opt falls relative to the universal constant, and whether the measurement budget can accommodate photon counting. Those two conditions, not the mode number n alone, determine which guarantee is available and at what cost.

Sources used in this answer