Proof of the hiding conjecture for Gaussian boson sampling with an arbitrary number of squeezed input modes

Laura Shou¹, Alexey V. Gorshkov^1,3, Victor Galitski¹, Sarah H. Miller² Address: ¹Joint Quantum Institute, Department of Physics, NIST/University of Maryland, College Park, MD 20742, USA Address: ²Applied Research Laboratory for Intelligence and Security, University of Maryland, College Park, Maryland 20742, USA Address: ³Joint Center for Quantum Information and Computer Science, NIST/University of Maryland, College Park, MD, 20742, USA

Abstract.

Gaussian boson sampling (GBS) is a sampling task proposed to demonstrate quantum advantage. We consider Gaussian boson sampling on $M$ optical modes, with $K$ equally squeezed input modes and $N$ observed photon counts. We complete the proof of the hiding conjecture for Gaussian boson sampling with an arbitrary number of squeezers $K$ , which is a part of the argument for classical hardness of GBS. In particular, we show that for any $K$ and $N=o(\sqrt{K})$ , the symmetric product $MK^{-1/2}U_{NK}U_{NK}^{T}$ , for $U_{NK}$ the top left $N\times K$ submatrix of an $M\times M$ Haar random unitary $U$ , is close in total variation distance to both an $N\times N$ symmetric complex Gaussian matrix $\mathbf{G}$ with independent entries, and the symmetric product $GG^{T}/\sqrt{K}$ for $G$ an $N\times K$ matrix of iid standard complex Gaussians. We show however that the density-based instance generating method of [AA13, Lemma 5.8] used to efficiently implement a hiding procedure fails for Gaussian boson sampling with $K=cM$ if $c<1/2$ . Instead we use approximate instance generating to implement the hiding for the usual classical hardness reduction.

1. Introduction

Gaussian boson sampling [HKS+17] is a sampling task which is expected to be hard for classical computers, but currently realizable in existing quantum experiments. In this task, one prepares an initial Gaussian state consisting of $M$ single-mode squeezed vacuum states with squeezing parameters $s_{i}\geq 0$ . For simplicity we take the first $K$ modes to have the same squeezing parameters $s_{i}=s>0$ , and the remaining $M-K$ modes to have the vacuum state. The initial state is then inserted into an $M$ -mode passive linear optical network described by an $M\times M$ linear optical unitary $U$ (Figure 1). After interfering in the linear optical network, the resulting output state is measured in the photon-number basis, producing a photon count outcome $\mathbf{n}\in\{0,1,2,\ldots\}^{M}$ , with total photon number $N=\sum_{i=1}^{M}\mathbf{n}_{i}$ . Let $I_{K}$ be the $M\times M$ diagonal matrix whose first $K$ diagonal entries are 1s, while the rest are 0s. Consider the $N\times N$ submatrix $(UI_{K}U^{T})_{\mathbf{n},\mathbf{n}}$ of the $M\times M$ matrix $UI_{K}U^{T}$ formed by taking the rows and columns corresponding to $1$ s in $\mathbf{n}$ . For collision-free outcomes $\mathbf{n}\in\{0,1\}^{M}$ , the probability of observing $\mathbf{n}$ is [HKS+17, KHS+19]

(1.1)

\displaystyle\mathbf{P}[\mathbf{n}]

\displaystyle=\frac{\tanh^{N}(s)}{\cosh^{K}(s)}|\operatorname{Haf}[(UI_{K}U^{T})_{\mathbf{n},\mathbf{n}}]|^{2},

where $\operatorname{Haf}[X]$ denotes the hafnian of a symmetric matrix $X$ ; letting $N=2n$ , the hafnian is $\operatorname{Haf}[X]=\sum_{\pi\in\mathcal{P}_{2}(2n)}\prod_{\{i,j\}\in\pi}X_{ij}$ , where $\mathcal{P}_{2}(2n)$ is the set of all pairings (i.e. perfect matchings) of $2n$ elements. Conditioned on observing total photon number $N$ , the collision-free condition for Haar random $U$ occurs with high probability if $N=o(M^{1/2})$ [AA13, DMV+22]. Since the average number of photons is $\langle N\rangle=K\sinh^{2}s$ , one takes the squeezing parameter $s$ to be small to ensure $\langle N\rangle=o(M^{1/2})$ .

Figure 1. Illustration of Gaussian boson sampling experiment.

K

single-mode squeezed states are inserted into an

M

-mode linear optical network described by an

M\times M

linear optical unitary

U

. The output state is measured in the photon number basis, producing a random sample

\mathbf{n}\in\{0,1,2,\ldots\}^{M}

of photon counts. For collision-free outcomes,

\mathbf{n}\in\{0,1\}^{M}

The argument that it is classically hard to generate samples $\mathbf{n}$ according to the distribution (1.1) is based on classical hardness of exactly computing permanents and hafnians of complex matrices [Val79]. However, since a realistic Gaussian boson sampler will always have some amount of noise and errors, one has to consider the task of approximate sampling from the distribution described by (1.1). The hardness of approximate average-case sampling then relies on two properties:

(1)

A complexity theoretic conjecture that approximating the hafnian of a random complex Gaussian-type matrix to a certain additive error is #P-hard in the average-case.
(2)

A hiding conjecture, which essentially states that one can “hide” a random complex Gaussian matrix $X$ as a submatrix of $UI_{K}U^{T}$ in total variation distance (Conjecture 1). This would allow one to use a GBS oracle running $U$ to approximate $|\operatorname{Haf}(X)|^{2}$ in $\text{BPP}^{\text{NP}}$ .

We focus on the hiding conjecture (2). Previously, only special cases of $K$ were proved to satisfy the hiding property, namely in the sparse squeezer regime $NK=o(M)$ [AA13, DMV+22, SMG26], and the case $K=M$ [SMG26]. Here we prove the hiding conjecture for all other $K$ , which proves the full range for the hiding conjecture for Gaussian boson sampling.

1.1. Preliminaries

To state the precise hiding conjecture and results, we first recall some definitions. The total variation distance (TVD) between two probability measures $\mu$ and $\nu$ on a measure space $(E,\mathcal{E})$ is

(1.2)

\displaystyle d_{\mathrm{TV}}(\mu,\nu)

\displaystyle\equiv\sup_{A\in\mathcal{E}}|\mu(A)-\nu(A)|.

It also has the characterization

(1.3)

\displaystyle d_{\mathrm{TV}}(\mu,\nu)

\displaystyle=\sup_{\|f\|_{\infty}\leq 1}\frac{1}{2}\int f\,(d\mu-d\nu).

If $\mu$ and $\nu$ have densities $p$ and $q$ with respect to a measure $d\lambda$ on $E$ , then also

(1.4)

\displaystyle d_{\mathrm{TV}}(\mu,\nu)

\displaystyle=\frac{1}{2}\int_{E}|p(x)-q(x)|\,d\lambda(x).

For random variables $Z$ and $W$ , we write $d_{\mathrm{TV}}(Z,W)$ to mean the total variation distance between their distributions $\mathcal{L}(Z)$ and $\mathcal{L}(W)$ .

For any coupling $(X,Y)$ of $(\mu,\nu)$ , i.e. random variables $X,Y$ such that $X\sim\mu$ and $Y\sim\nu$ , there is the TVD coupling bound $d_{\mathrm{TV}}(\mu,\nu)\leq\mathbf{P}[X\neq Y]$ .

We will use the following probability distributions.

•

Let $\mathcal{CN}(0,\sigma^{2})$ denote the complex Gaussian distribution whose real and imaginary parts are independent Gaussians with mean $0$ and variance $\sigma^{2}/2$ .
•

Let $\mathcal{G}_{N}^{\mathrm{sym}}$ denote the ensemble of $N\times N$ symmetric random matrices $\mathbf{G}$ with $\mathcal{CN}(0,2)$ diagonal entries and $\mathcal{CN}(0,1)$ off-diagonal entries, with all entries independent modulo the symmetry requirement. We will use bold font $\mathbf{G}$ only for $\mathbf{G}\sim\mathcal{G}_{N}^{\mathrm{sym}}$ .
•

Let $\mathcal{G}\mathcal{G}^{T}_{NK}$ be the ensemble of $N\times N$ matrices $GG^{T}/\sqrt{K}$ where $G$ is an $N\times K$ matrix of iid standard complex Gaussian entries. The normalization is chosen so that $GG^{T}/\sqrt{K}$ has a nondegenerate limiting distribution as $K\to\infty$ . Note that since $G$ is complex, $GG^{T}$ is not Wishart.
•

Let $\mathcal{H}_{M}$ denote the Haar measure on the unitary group $\mathrm{U}(M)$ .

The hiding conjecture can be stated as follows [AA13, DMV+22, SMG26]:

Conjecture 1 (hiding in Gaussian boson sampling).

Let $U_{NK}$ be the top left $N\times K$ submatrix of an $M\times M$ Haar random unitary matrix, and let $Z=Z_{N,K}$ be a matrix with distribution given by either $\mathcal{G}_{N}^{\mathrm{sym}}$ or $\mathcal{G}\mathcal{G}^{T}_{NK}$ . Then for $N\leq K\leq M$ , there exist polynomials $p,r$ such that for any $\delta>0$ and $M\geq p(N)/r(\delta)$ , at least one of the choices of $Z$ satisfies

(1.5)

\displaystyle d_{\mathrm{TV}}(MK^{-1/2}U_{NK}U_{NK}^{T},Z)=O(\delta),

where $d_{\mathrm{TV}}(\cdot,\cdot)$ denotes the total variation distance as defined in (1.2).

The sparse squeezer case $N=K=o(M^{1/5})$ was observed [HKS+17, DMV+22] to follow from hiding in Fock boson sampling [AA13] with $Z\sim\mathcal{G}\mathcal{G}^{T}_{NK}$ . However, the non-sparse regime where $K$ is large is of most interest experimentally [ZWD+20, ZDQ+21, MLA+22, DGL+23, LSD+26]. Additionally, a large number of squeezers $K$ is favorable for the anticoncentration results of [EID+25b, EID+25a], which provide evidence for sampling hardness. In the large $K\propto M$ regime, only the case $K=M$ was proved, with $Z\sim\mathcal{G}_{N}^{\mathrm{sym}}$ , using that in this case the matrix $U_{NK}U_{NK}^{T}$ is a submatrix of a COE (Circular Orthogonal Ensemble) random matrix [SMG26].

1.2. Main results

In this paper, we prove the hiding conjecture for Gaussian boson sampling for any number of input squeezed modes $K$ , where $N\leq K\leq M$ . This fully resolves Conjecture 1. Additionally, we prove TVD closeness of the distributions $\mathcal{G}_{N}^{\mathrm{sym}}$ and $\mathcal{G}\mathcal{G}^{T}_{NK}$ in the regime $N=o(\sqrt{K})$ , so that Conjecture 1 holds with either distribution in this regime.

Recall in the context of Gaussian boson sampling, $M$ is the number of modes, $K$ is the number of squeezed modes, and $N$ is the number of detected photons. While $N$ is even for Gaussian boson sampling, we do not require $N$ to be even in Theorems 1.1 and 1.2 below.

Theorem 1.1 (hiding in Gaussian boson sampling).

Let $N\leq K\leq M$ , let $U_{NK}$ be the top left $N\times K$ submatrix of an $M\times M$ Haar random unitary matrix $U$ , and let $\mathbf{G}\sim\mathcal{G}_{N}^{\mathrm{sym}}$ . Then as $K\to\infty$ ,

(1.6)

\displaystyle d_{\mathrm{TV}}(MK^{-1/2}U_{NK}U_{NK}^{T},\mathbf{G})=O\left(\frac{N}{\sqrt{K}}\right),

which is $o(1)$ if $N=o(\sqrt{K})$ .

Remark 1.1.

(1)
Theorem 1.1 combined with the sparse case result [SMG26, Theorem 1.5] for $NK=o(M)$ gives the full range of Conjecture 1 by choosing polynomials $p,r$ appropriately. For example:
- •
  
  For $K\propto M$ , one can take $Z\sim\mathcal{G}_{N}^{\mathrm{sym}}$ , $p(N)=N^{2}$ , and $r(\delta)=c\delta^{2}$ for sufficiently small $c$ , by Theorem 1.1.
- •
  
  For $K\leq M^{1/2}$ , one can take $Z\sim\mathcal{G}\mathcal{G}^{T}_{NK}$ , $p(N)=N^{2+\epsilon}$ , and $r(\delta)=c\delta^{4}$ by the sparse result, rewritten as (3.1).
- •
  
  For any $K\geq M^{1/2}$ , one can take $Z\sim\mathcal{G}_{N}^{\mathrm{sym}}$ , $p(N)=N^{4}$ , and $r(\delta)=c\delta^{4}$ , again by Theorem 1.1.
In fact, using Theorem 1.2 below, we also obtain Conjecture 1 with a fixed distribution $Z=GG^{T}/\sqrt{K}\sim\mathcal{G}\mathcal{G}^{T}_{NK}$ and fixed polynomials $p,r$ , over any $N\leq K\leq M$ : For $p(N)=N^{4}$ , $r(\delta)=c\delta^{4}$ , and any $M\geq p(N)/r(\delta)$ ,

(1.7) $\displaystyle d_{\mathrm{TV}}(MK^{-1/2}U_{NK}U_{NK}^{T},GG^{T}/\sqrt{K})=O(\delta).$
(2)

In the main regime of interest $K\propto M$ , the condition $N=o(\sqrt{K})$ becomes $N=o(\sqrt{M})$ , which is the conjectured maximal submatrix size allowed [DMV+22]. As we will show in Theorem 1.2, in the regime $N=o(\sqrt{K})$ , TVD closeness in (1.5) actually holds with either distribution $Z\sim\mathcal{G}_{N}^{\mathrm{sym}}$ or $\mathcal{G}\mathcal{G}^{T}_{NK}$ .
(3)

The proof of Theorem 1.1 relies on relating the $K<M$ case to the $K=M$ case. One could alternatively write down the density formula for $U_{NK}U_{NK}^{T}$ in terms of a matrix integral and use quantitative Laplace approximation to bound the TVD. However, due to error bounds for high-dimensional Laplace’s method, we do not expect this approach, at least with standard estimates, to obtain any sharper results.

Let $G$ be an $N\times K$ matrix of iid standard complex Gaussians. Combining Theorem 1.1 with TVD closeness of $MU_{NK}U_{NK}^{T}$ to $GG^{T}$ in a “sparse squeezer” regime with $N\leq K$ and $NK=o(M)$ [SMG26], we will obtain

Theorem 1.2 ( $GG^{T}$ vs $\mathcal{G}_{N}^{\mathrm{sym}}$ ).

Let $G$ be an $N\times K$ matrix of iid standard complex Gaussians, and let $\mathbf{G}\sim\mathcal{G}_{N}^{\mathrm{sym}}$ . Then as $K\to\infty$ ,

(1.8)

\displaystyle d_{\mathrm{TV}}(GG^{T}/\sqrt{K},\mathbf{G})=O\left(\frac{N}{\sqrt{K}}\right),

which is $o(1)$ if $N=o(\sqrt{K})$ .

This implies that when $N=o(\sqrt{K})$ , it does not matter whether we use $GG^{T}$ or $\mathbf{G}\sim\mathcal{G}_{N}^{\mathrm{sym}}$ as the target Gaussian-type distribution in the hiding statement. We note that the matrices $\mathbf{G}$ are much nicer to work with due to their independent entries. Moreover, they are much closer to the random Gaussian matrices used in the hardness reduction for Fock boson sampling, which suggests techniques and results for hardness of Fock boson sampling will carry over more easily to $\mathbf{G}$ than to $GG^{T}$ .

Remark 1.2.

The $N=o(\sqrt{K})$ condition differs from the $\Theta(K^{1/3})$ transition for GOE (Gaussian Orthogonal Ensemble) behavior of Wishart matrices with $K$ degrees of freedom [BDER16, JL15, RR19]. It remains to locate the precise location of the transition in this case.

The TVD hiding property Theorem 1.1 demonstrates it is possible to “hide” a Gaussian matrix $\mathbf{G}\sim\mathcal{G}_{N}^{\mathrm{sym}}$ as a random instance of $MK^{-1/2}(UI_{K}U^{T})_{S,S}$ for $U$ a Haar random unitary matrix and $S$ a random size $N$ subset of $\{1,\ldots,M\}$ , up to small TVD error. However, for the usual classical hardness argument, given $\mathbf{G}$ , we need to generate an instance of such a random $U$ and $S^{*}$ efficiently. The instance-generating procedure for Fock boson sampling uses a rejection sampling method based on the density bound $\tilde{f}(Z)\leq(1+o(1))\tilde{g}(Z)$ , where $\tilde{f}$ is the unitary submatrix density and $\tilde{g}$ the Gaussian density [AA13, Lemma 5.7]. However, we show that such a pointwise density bound cannot hold for Gaussian boson sampling with $K=\alpha M$ , $0<\alpha<1/2$ .

Proposition 1.3 (no density ratio bound).

Let $K\geq 4$ , $K\geq N$ , and $N+K\leq M$ . Denote by $U_{NK}$ the top left $N\times K$ submatrix of an $M\times M$ Haar random unitary matrix $U$ , and let $\mathbf{G}\sim\mathcal{G}_{N}^{\mathrm{sym}}$ . Let $f$ be the density function of $MK^{-1/2}U_{NK}U_{NK}^{T}$ (if it exists) and $g$ be the density function of $\mathbf{G}$ over the space of $N\times N$ complex symmetric matrices¹¹ 1 We view this as the space spanned by the upper triangular matrix elements. Then for $K=\alpha M$ , $0<\alpha<1/2$ , and any $N\leq K$ , either the density function doesn’t exist, or

(1.9)

\displaystyle\operatorname{ess\,sup}_{Z}\frac{f(Z)}{g(Z)}\geq\Omega_{\alpha}(e^{c_{\alpha}M}),

for a constant $c_{\alpha}>0$ and where $\Omega_{\alpha}$ indicates the implicit constant may also depend on $\alpha$ .

In particular, consider a sequence of $K=K(M),N=N(M)$ with $\alpha=\alpha(M)=K(M)/M$ bounded in $[\eta,1/2-\eta]$ for some $\eta>0$ , and let $f_{M,K,N},g_{M,K,N}$ denote the density functions as above if they exist. Then

(1.10)

\displaystyle\liminf_{M\to\infty}\operatorname{ess\,sup}_{Z}\frac{f_{M,K,N}(Z)}{g_{M,K,N}(Z)}=+\infty.

Remark 1.3.

This result is perhaps counterintuitive, since intuitively, smaller $K\propto M$ should not make the submatrix $U_{NK}U_{NK}^{T}$ look less Gaussian than the $K=M$ case, which sees the full orthonormality requirements of $U$ . In the $K=M$ case, the bound $f(Z)\leq(1+o(1))g(Z)$ holds [SMG26] for $N=o(K^{1/3})$ . However, the pointwise density bound (1.9) captures rare tail behavior, which is not necessarily captured by TVD, and which plausibly can differ for smaller $K$ where $MK^{-1/2}U_{NK}U_{NK}^{T}$ may have different tail behavior.

Proposition 1.3 means we cannot use the density-based rejection sampling method of [AA13, Lemmas 5.7, 5.8], which requires $f(Z)\leq(1+\delta)g(Z)$ , to generate instances of unitaries $U$ with $\mathbf{G}$ hidden as an instance of $MK^{-1/2}(UI_{K}U^{T})_{S^{*},S^{*}}$ when $K<M/2$ . One could perhaps try to prove the required density bound for $K>M/2$ (we know it at least holds for $K=M$ and $N=o(K^{1/3})$ ), or try to truncate the distribution $\mathbf{G}$ to remove tail behavior, but both of these approaches are model-specific and also likely involve working with complicated density functions. Instead, we adapt the argument of [AA13, §5.2] to generate approximately Haar instances of unitaries $U$ and approximately uniform locations $S$ to hide $\mathbf{G}$ , using postselection with an NP oracle. This will be enough to complete the hardness reduction, under a finite-precision implementation assumption, giving Theorem 1.4 below. To state the theorem, first define

Problem 1 ( $|\mathrm{GHE}|_{\pm}^{2}$ ).

Let $N\in 2\mathbb{N}$ . Given as input a matrix $X\sim\mathcal{G}_{N}^{\mathrm{sym}}$ , together with error bounds $\varepsilon,\delta>0$ , estimate $|\operatorname{Haf}(X)|^{2}$ to within additive error $\pm\varepsilon\cdot\mathbf{E}|\operatorname{Haf}(X)|^{2}$ with probability at least $1-\delta$ over $X$ and the algorithm’s randomness in $\operatorname{poly}(N,1/\varepsilon,1/\delta)$ time.

For $X\sim\mathcal{G}_{N}^{\mathrm{sym}}$ , one can use independence of entries above the diagonal to quickly calculate that

(1.11)

\displaystyle\mathbf{E}|\operatorname{Haf}(X)|^{2}

\displaystyle=(N-1)!!=\frac{N!}{(N/2)!2^{N/2}}\sim\frac{\sqrt{2}N^{N/2}}{e^{N/2}}.

As is usual [AA13, §2], it will be understood that all entries of $X$ are rounded to polynomially many bits of precision; here we allow $\operatorname{poly}(M,1/\delta,1/\varepsilon)$ bits of precision, for a sufficiently large fixed polynomial. To formalize this in our case, we will make an assumption on finite-precision implementation, stated later precisely as Assumption 1 in Section 4. The assumption essentially says that we can approximate Haar random $U$ by finite-precision descriptions $v=v_{\xi}$ , and also $(UI_{K}U^{T})_{S,S}$ by an efficiently computable finite-precision matrix $\hat{Z}$ , using high enough numerical precision. We expect this assumption is true and that it can be proved by careful finite-precision accounting and Gram-Schmidt or QR factorization.

Under this finite-precision implementation assumption, we prove the analogue of the main boson sampling result of [AA13, Theorem 1.3], for Gaussian boson sampling with essentially arbitrary number of squeezed input modes $K$ . We refer to [com] for definitions of the standard complexity classes NP, BPP (bounded-error probabilistic polynomial-time), and FBPP (the function/search analogue of BPP, which searches for a witness to a relation in probabilistic polynomial time). The notation $\text{FBPP}^{\text{NP}}$ means FBPP with access to an NP oracle.

Theorem 1.4 (main hardness result).

Let the probability distribution $\mathcal{D}_{A}$ be the output of a Gaussian boson sampling experiment $A=A(v,K,s)$ , for $v$ a finite-precision description of a linear optical unitary $U=U(v)$ , and $K,s$ input squeezing parameters. Suppose there exists a classical algorithm $C$ which is able to approximately sample from limited instances of $A$ , say for each $M\in\mathbb{N}$ , only those with $K=K(M)$ equally squeezed input modes for some given²² 2 We assume the sequence $K(M)$ is efficiently computable. sequence $K(M)$ with $M=O(\operatorname{poly}(K(M)))$ , and all small squeezing parameters $s=O(K^{-1/4})$ . More precisely, the algorithm $C$ takes as input a description of such $A$ as well as an error bound $\varepsilon$ , and samples from a probability distribution $\mathcal{D}^{\prime}_{A}$ such that $\|\mathcal{D}^{\prime}_{A}-\mathcal{D}_{A}\|_{\mathrm{TV}}\leq\varepsilon$ in $\operatorname{poly}(|A|,1/\varepsilon)$ time, where $|A|$ denotes the length of the description³³ 3 As in [AA13], and as discussed more in Section 4, it will be understood that all entries in $A$ are rounded to $\operatorname{poly}(M,1/\delta,1/\varepsilon)$ bits of precision. of $A$ . Then under the finite-precision implementation Assumption 1, the problem $|\mathrm{GHE}|_{\pm}^{2}$ is solvable in $\text{FBPP}^{\text{NP}}$ . In other words, if we treat $C$ as a black box, then $|\mathrm{GHE}|_{\pm}^{2}\in\text{FBPP}^{\text{NP}^{C}}$ .

Conjecture 2.

$|\mathrm{GHE}|_{\pm}^{2}$ is #P-hard, in the sense that if $\mathcal{O}$ is any oracle that solves $|\mathrm{GHE}|_{\pm}^{2}$ , then $\text{P}^{\#\text{P}}\subseteq\text{BPP}^{\mathcal{O}}$ .

Conjecture 2 is analogous to the conjecture for additive approximation of squared permanents of complex matrices, $|\mathrm{GPE}|_{\pm}^{2}$ , being #P-hard [AA13]. Note that in the formulation of $|\mathrm{GHE}|_{\pm}^{2}$ here as well as in $|\mathrm{GPE}|_{\pm}^{2}$ , the additive error size is $\varepsilon$ times the average value of the quantity to estimate ( $|\operatorname{Haf}(X)|^{2}$ here, $|\operatorname{Per}(X)|^{2}$ for $|\mathrm{GPE}|_{\pm}^{2}$ ). This makes Problem 1 the natural hafnian analogue of $|\mathrm{GPE}|_{\pm}^{2}$ from [AA13] (previous hafnian hardness problems in the context of GBS used the more complicated matrix $XX^{T}$ for $X$ an $N\times K$ matrix of iid complex Gaussians, and also did not express the error bound in terms of the hafnian moments). If Conjecture 2 holds, then the existence of such a classical algorithm $C$ in Theorem 1.4 implies collapse of the polynomial hierarchy by Toda’s theorem [Tod91].

1.3. Outline

The rest of the paper is organized as follows. Additionally, we note that generic constants $C,c$ may change from line to line throughout the paper.

•

In Section 2, we prove a simpler version of the hiding property Theorem 1.1, though which holds only for $N=o(K^{1/3})$ in general. This proof is less technical, and is enough to perform the subsequent hardness reduction. The extension to $N=o(\sqrt{K})$ is proved in Appendix A.
•

In Section 3, we prove Theorem 1.2 on $\mathcal{G}_{N}^{\mathrm{sym}}$ vs $\mathcal{G}\mathcal{G}^{T}_{NK}$ .
•

In Section 4, we prove the hardness reduction Theorem 1.4.
•

In Section 5, we prove Proposition 1.3 on the lack of a density ratio bound for $K=\alpha M$ , $0<\alpha<1/2$ .
•

In Appendix A, we prove the full Theorem 1.1 up to submatrix size $N=o(\sqrt{K})$ .

2. Proof of $N=o(K^{1/3})$ TVD hiding

In this section we prove a weaker version of Theorem 1.1; namely,

Proposition 2.1 (hiding for $N=o(K^{1/3})$ ).

(2.1)

\displaystyle d_{\mathrm{TV}}(MK^{-1/2}U_{NK}U_{NK}^{T},\mathbf{G})

\displaystyle=O\left(\sqrt{\frac{N^{3}}{K}}\right),

which is $o(1)$ if $N=o(K^{1/3})$ .

This is strictly weaker than the bound and $N=o(K^{1/2})$ allowed in Theorem 1.1, but the proof of Proposition 2.1 will be simpler, and this statement is enough to go through the usual hardness reduction with some minor parameter adjustments. (Recall, the original boson sampling hardness reduction of [AA13] had $N=o(M^{1/5})$ ; proving TVD closeness with any inverse polynomial power is sufficient.) The key idea of the proof of Proposition 2.1 or Theorem 1.1 is Lemma 2.2, which is presented in this section. The proof of the stronger bound in Theorem 1.1 differs only in a later technical estimate, which we give in Appendix A.

It suffices to prove (2.1) for $N^{3}=O(K)$ ; otherwise the bound is trivial. Let $\mathcal{V}$ be the span of the last $M-K$ columns of $U$ . The main idea is the following: The first $K$ columns of $U$ form an orthonormal basis for $\mathcal{V}^{\perp}$ , and conditioned on $\mathcal{V}$ they form a Haar random basis for $\mathcal{V}^{\perp}$ [Mec19, §1.2]. Since $\mathcal{V}^{\perp}$ has dimension $K$ , this can be used to relate the distribution of $U_{NK}U_{NK}^{T}$ to the distribution of a matrix transformation of $W_{NK}W_{NK}^{T}$ , where $W$ is a $K\times K$ Haar random unitary [Eq. (2.2)]. The matrix $W$ can be thought of as acting as the source of randomness for the random basis of $\mathcal{V}^{\perp}$ . The distribution of $W_{NK}W_{NK}^{T}$ is the setting considered in [SMG26], and so it will be possible to obtain Theorem 1.1 for $K<M$ by reducing to the hiding case with $K=M$ proved there.

We start with the first part, on relating the distribution of $U_{NK}U_{NK}^{T}$ to one involving $W_{NK}W_{NK}^{T}$ for $W$ a $K\times K$ Haar random unitary.

Lemma 2.2.

Let $U$ be an $M\times M$ Haar random matrix, and let $U_{NK}$ be its top left $N\times K$ submatrix with $N\leq K$ . Denote by $V$ the $M\times(M-K)$ matrix consisting of the last $M-K$ columns of $U$ , and define the $N\times N$ matrix $P_{0}:=I_{N}-E_{N}VV^{\dagger}E_{N}^{\dagger}$ , for $E_{N}=(I_{N}\;0_{M-N})$ . Then

(2.2)

\displaystyle U_{NK}\overset{d}{=}P_{0}^{1/2}W_{NK},\quad\text{ and }\quad U_{NK}U_{NK}^{T}\overset{d}{=}P_{0}^{1/2}W_{NK}W_{NK}^{T}(P_{0}^{1/2})^{T},

for $W$ an independent $K\times K$ Haar random unitary and $W_{NK}$ its top $N\times K$ submatrix.

Proof.

Let $\mathcal{V}$ be the span of the last $M-K$ columns of $U$ . Let $A_{0}$ be an $M\times K$ matrix whose columns form an orthonormal basis for $\mathcal{V}^{\perp}$ , and let $W$ be an independent $K\times K$ Haar random unitary. Conditioned on $\mathcal{V}$ , we have

(2.3)

\displaystyle U_{NK}\overset{d}{=}E_{N}A_{0}W,

for $E_{N}=(I_{N}\;0_{M-N})$ . Since the distribution of $W$ is invariant under unitary rotations, we would like to pass the $E_{N}$ through to $W$ (so that we can consider $W_{NK}$ ), but $A_{0}$ is rectangle-shaped so we have to do some manipulation.

Note that $A_{0}A_{0}^{\dagger}=I_{M}-VV^{\dagger}$ the projection onto $\mathcal{V}^{\perp}$ , and so $(E_{N}A_{0})(E_{N}A_{0})^{\dagger}=I_{N}-E_{N}VV^{\dagger}E_{N}^{\dagger}=P_{0}$ . We can do polar decomposition (or SVD) to obtain $E_{N}A_{0}=P_{0}^{1/2}R$ , for $R$ an $N\times K$ semiunitary matrix;⁴⁴ 4 In general, $R$ need not be unique. However, in this case one can show $P_{0}$ is invertible almost surely (a.s.) since $N\leq K$ , so $R$ is unique a.s. To see this, let $A$ be the $M\times K$ matrix consisting of the first $K$ columns of $U$ , so $P_{0}=(E_{N}A)(E_{N}A)^{\dagger}$ . For $Z$ an $M\times K$ matrix of iid standard complex Gaussians, we have $E_{N}A\overset{d}{=}E_{N}Z(Z^{\dagger}Z)^{-1/2}$ [Mec19, §1.2]. Note that $Z$ has rank $K$ a.s. (the probability of the next column being in the span of the previous columns is 0), so $\operatorname{rank}Z^{\dagger}Z=\operatorname{rank}Z=K$ and $Z^{\dagger}Z$ is invertible a.s. The matrix $E_{N}Z$ is an $N\times K$ matrix of iid standard complex Gaussians, so has rank $N$ a.s. Thus $\operatorname{rank}[(E_{N}A)(E_{N}A)^{\dagger}]=\operatorname{rank}(E_{N}A)=N$ a.s., so $P_{0}$ is invertible a.s. i.e. $RR^{\dagger}=I_{N}$ and the rows of $R$ are orthonormal. Then we can extend the rows of $R$ to a full orthonormal basis of $\mathbb{C}^{K}$ , which will be given by the rows of a $K\times K$ unitary $\xi$ (independent of $W$ ), with $R=(I_{N}\;0_{K-N})\xi$ . Since $W$ is invariant under unitary multiplication, and $\xi$ is independent of $W$ , conditioned on $\mathcal{V}$ we get

	$\displaystyle U_{NK}$	$\displaystyle\overset{d}{=}E_{N}A_{0}W=P_{0}^{1/2}(I_{N}\;0_{K-N})\xi W$
(2.4)			$\displaystyle\overset{d}{=}P_{0}^{1/2}(I_{N}\;0_{K-N})W=P_{0}^{1/2}W_{NK},$

where $W_{NK}$ is the top $N\times K$ submatrix of the $K\times K$ unitary $W$ . Thus we obtain (2.2). ∎

Proof of Proposition 2.1.

By Lemma 2.2, the distribution of $U_{NK}U_{NK}^{T}$ is the same as the distribution of $P_{0}^{1/2}W_{NK}W_{NK}^{T}(P_{0}^{1/2})^{T}$ , where $P_{0}$ is defined as in the lemma and $W$ is an independent $K\times K$ Haar random unitary matrix. By [SMG26], since $W_{NK}$ is the $N\times K$ submatrix for a $K\times K$ Haar unitary matrix, then for $\mathbf{G}\sim\mathcal{G}_{N}^{\mathrm{sym}}$ and $N=O(\sqrt{K})$ ,

(2.5)

\displaystyle d_{\mathrm{TV}}(\sqrt{K}W_{NK}W_{NK}^{T},\mathbf{G})\leq O(N/\sqrt{K}).

Letting $Q:=(MK^{-1}P_{0})^{1/2}$ , which is independent of $\mathbf{G}$ and $W$ , we can then write

	$\displaystyle d_{\mathrm{TV}}(MK^{-1/2}U_{NK}U_{NK}^{T},\mathbf{G})$	$\displaystyle\leq d_{\mathrm{TV}}(Q\sqrt{K}W_{NK}W_{NK}^{T}Q^{T},Q\mathbf{G}Q^{T})+d_{\mathrm{TV}}(Q\mathbf{G}Q^{T},\mathbf{G})$
(2.6)			$\displaystyle\leq O(N/\sqrt{K})+d_{\mathrm{TV}}(Q\mathbf{G}Q^{T},\mathbf{G}).$

So to prove the theorem it suffices to bound $d_{\mathrm{TV}}(Q\mathbf{G}Q^{T},\mathbf{G})$ . To do this, we will show that $Q$ is typically close to the identity, which will make the TVD small.

For fixed $Q=q$ (e.g. conditioned on $V$ ), both of the involved distributions $q\mathbf{G}q^{T}$ and $\mathbf{G}$ have explicit Gaussian densities, so we can directly estimate the TVD. For fixed invertible $q$ , the density of $q\mathbf{G}q^{T}$ over the space of $N\times N$ symmetric complex matrices is calculated by change of variables⁵⁵ 5 The map $X\mapsto qXq^{T}$ is linear, and can be expressed as $(q\otimes q)\operatorname{vec}(X)\equiv\operatorname{vec}(qXq^{T})$ , where $\operatorname{vec}(X)$ is the vectorization of $X$ formed by stacking columns of $X$ . To calculate the Jacobian of the map $X\mapsto qXq^{T}$ , if $q$ is diagonalizable with eigenpairs $\{(\lambda_{i},|u_{i}\rangle)\}_{i=1}^{N}$ , take the eigenbasis $u_{i}u_{j}^{T}+u_{j}u_{i}^{T}=|u_{i}\rangle\langle\bar{u}_{j}|+|u_{j}\rangle\langle\bar{u}_{i}|$ for $i\leq j$ of the transformation $X\mapsto qXq^{T}$ , which gives (complex) Jacobian $\prod_{i\leq j}\lambda_{i}\lambda_{j}=(\det q)^{N+1}$ and (real) Jacobian $|\det q|^{2(N+1)}$ . from $\mathbf{G}$ to be

(2.7)

\displaystyle f_{q}(Z)=c_{N}|\det q|^{-2(N+1)}e^{-\frac{1}{2}\operatorname{Tr}\bar{q}^{-1}Z^{\dagger}(q^{-1})^{\dagger}q^{-1}Z(q^{-1})^{T}},

where $c_{N}=2^{-N}\pi^{-N(N+1)/2}$ is the normalization constant for the density $c_{N}e^{-\frac{1}{2}\operatorname{Tr}Z^{\dagger}Z}$ of $\mathbf{G}$ .

Total variation distance can be bounded using the Kullback–Leibler (KL) divergence, or relative entropy, via Pinsker’s inequality,

(2.8)

\displaystyle d_{\mathrm{TV}}(X,Y)

\displaystyle\leq\sqrt{\frac{1}{2}D_{\mathrm{KL}}(X||Y)},

where the KL divergence is $D_{\mathrm{KL}}(X||Y):=\int f(x)\log\frac{f(x)}{g(x)}\,dx$ for $f$ and $g$ the respective density functions for $X$ and $Y$ .

For fixed invertible $q$ , letting $s=q^{\dagger}q$ , the KL divergence for $q\mathbf{G}q^{T}$ and $\mathbf{G}$ is

	$\displaystyle D_{\mathrm{KL}}(q\mathbf{G}q^{T}\|\|\mathbf{G})$	$\displaystyle=\mathbf{E}_{Z\sim q\mathbf{G}q^{T}}[\log\|\det q\|^{-2(N+1)}-\frac{1}{2}\operatorname{Tr}\bar{q}^{-1}Z^{\dagger}(q^{-1})^{\dagger}q^{-1}Z(q^{-1})^{T}+\frac{1}{2}\operatorname{Tr}Z^{\dagger}Z]$
(2.9)			$\displaystyle=-(N+1)\log\det s-\frac{1}{2}[N^{2}+N]+\frac{1}{2}[(\operatorname{Tr}s)^{2}+\operatorname{Tr}(s^{2})],$

using e.g. $\mathbf{E}_{Y\sim\mathcal{G}_{N}^{\mathrm{sym}}}[\operatorname{Tr}(Y^{\dagger}sYs^{T})]=(\operatorname{Tr}s)^{2}+\operatorname{Tr}(s^{2})$ by direct expansion of the trace. Write $s=q^{\dagger}q=I_{N}+\delta$ . (We will later show that, for random $Q$ , $\delta$ is small with high probability, so $s$ is typically a small perturbation of the identity.) Expand

	$\displaystyle D_{\mathrm{KL}}(q\mathbf{G}q^{T}\|\|\mathbf{G})$	$\displaystyle=-(N+1)\log\det(I_{N}+\delta)-\frac{1}{2}[N^{2}+N]+\frac{1}{2}[N^{2}+2N\operatorname{Tr}\delta+(\operatorname{Tr}\delta)^{2}+N+2\operatorname{Tr}\delta+\operatorname{Tr}(\delta^{2})]$
		$\displaystyle=(N+1)[\operatorname{Tr}\delta-\operatorname{Tr}\log(I_{N}+\delta)]+\frac{1}{2}(\operatorname{Tr}\delta)^{2}+\frac{1}{2}\operatorname{Tr}(\delta^{2})$
(2.10)			$\displaystyle\leq(N+3/2)\operatorname{Tr}(\delta^{2})+\frac{1}{2}(\operatorname{Tr}\delta)^{2},\quad\text{ if }\\|\delta\\|\leq 1/2,$

using $x-\log(1+x)\leq x^{2}$ for $-1/2\leq x$ .

Now we return to random $P_{0}$ and $Q=(MK^{-1}P_{0})^{1/2}$ . Let $E_{N}V=:v$ be the top right $N\times(M-K)$ submatrix of $U$ . We will show that $S=Q^{\dagger}Q=MK^{-1}(I_{N}-vv^{\dagger})$ is typically a small perturbation of the identity. Letting $L:=M-K$ for notational convenience, then

\displaystyle S=\frac{M}{K}I_{N}-\frac{M}{K}vv^{\dagger}

\displaystyle=I_{N}+\frac{L}{K}I_{N}-\frac{M}{K}vv^{\dagger}.

Then $\delta=\frac{L}{K}I_{N}-\frac{M}{K}vv^{\dagger}$ , and so $\operatorname{Tr}\delta=\frac{LN}{K}-\frac{M}{K}\operatorname{Tr}(vv^{\dagger})$ , and $\delta^{2}=\frac{L^{2}}{K^{2}}I_{N}-\frac{2LM}{K^{2}}vv^{\dagger}+\frac{M^{2}}{K^{2}}vv^{\dagger}vv^{\dagger}$ . Taking the expectation over the Haar random unitary $U$ , we see that $\mathbf{E}\operatorname{Tr}\delta=0$ since $\mathbf{E}\operatorname{Tr}(vv^{\dagger})=\frac{LN}{M}$ . Also, using Weingarten calculus, see e.g. [CS06, Col03],

	$\displaystyle\mathbf{E}(\operatorname{Tr}(vv^{\dagger}))^{2}=\sum_{i,k=1}^{N}\sum_{j,\ell=1}^{L}\mathbf{E}\|v_{ij}\|^{2}\|v_{k\ell}\|^{2}$	$\displaystyle=\sum_{i,k=1}^{N}\sum_{j,\ell=1}^{L}\frac{1}{M^{2}-1}[1+\delta_{ik}\delta_{j\ell}]-\frac{1}{M(M^{2}-1)}[\delta_{j\ell}+\delta_{ik}]$
(2.11)			$\displaystyle=\frac{N^{2}L^{2}+NL}{M^{2}-1}-\frac{N^{2}L+NL^{2}}{M(M^{2}-1)},$

and

	$\displaystyle\mathbf{E}\operatorname{Tr}(vv^{\dagger}vv^{\dagger})=\sum_{x,y=1}^{N}\sum_{i,j=1}^{L}\mathbf{E}[v_{xi}v_{yj}\bar{v}_{xj}\bar{v}_{yi}]$	$\displaystyle=\sum_{x,y=1}^{N}\sum_{i,j=1}^{L}\frac{1}{M^{2}-1}[\delta_{ij}+\delta_{xy}]-\frac{1}{M(M^{2}-1)}[1+\delta_{xy}\delta_{ij}]$
(2.12)			$\displaystyle=\frac{N^{2}L+NL^{2}}{M^{2}-1}-\frac{N^{2}L^{2}+NL}{M(M^{2}-1)}.$

Then

	$\displaystyle\mathbf{E}(\operatorname{Tr}\delta)^{2}$	$\displaystyle=\frac{L^{2}N^{2}}{K^{2}}-\frac{2LMN}{K^{2}}\mathbf{E}\operatorname{Tr}(vv^{\dagger})+\frac{M^{2}}{K^{2}}\mathbf{E}(\operatorname{Tr}(vv^{\dagger}))^{2}$
(2.13)			$\displaystyle=\frac{LN(M-N)}{K(M^{2}-1)}=O\left(\frac{LN}{KM}\right),$
	$\displaystyle\mathbf{E}\operatorname{Tr}(\delta^{2})$	$\displaystyle=\frac{L^{2}N}{K^{2}}-\frac{2LM}{K^{2}}\mathbf{E}\operatorname{Tr}(vv^{\dagger})+\frac{M^{2}}{K^{2}}\mathbf{E}\operatorname{Tr}(vv^{\dagger}vv^{\dagger})$
(2.14)			$\displaystyle=\frac{LN(MN-1)}{K(M^{2}-1)}=O\left(\frac{LN^{2}}{KM}\right).$

This also implies

(2.15)

\displaystyle\mathbf{P}[\|\delta\|>1/2]

\displaystyle\leq 4\mathbf{E}\operatorname{Tr}(\delta^{2})=O\left(\frac{LN^{2}}{KM}\right).

Returning to random $Q$ , in order to invoke (2.10), we need to ensure $\|\delta\|\leq 1/2$ , which occurs with probability at least $1-O(N^{2}/K)$ by (2.15). For handling the event $\|\delta\|>1/2$ , we want to go back to TVD since it is bounded, while KL-divergence need not be. Let $\delta(q):=q^{\dagger}q-I_{N}$ . Define a random variable $\tilde{Q}$ via

\displaystyle\tilde{Q}=\begin{cases}Q,&\|\delta(Q)\|\leq 1/2\\ I,&\text{otherwise}\end{cases};

then $d_{\mathrm{TV}}(Q,\tilde{Q})\leq\mathbf{P}[Q\neq\tilde{Q}]=\mathbf{P}[\|\delta(Q)\|>1/2]$ . Since $\tilde{Q}^{\dagger}\tilde{Q}=I_{N}+\delta(\tilde{Q})$ and $\delta(\tilde{Q})=\delta(Q)\mathbf{1}_{\|\delta(Q)\|\leq 1/2}$ , then

(2.16)

\displaystyle\begin{aligned} \mathbf{E}_{\tilde{Q}}[(\operatorname{Tr}\delta)^{2}]&\leq\mathbf{E}_{Q}[(\operatorname{Tr}\delta)^{2}]=O\left(\frac{LN}{KM}\right),\\ \mathbf{E}_{\tilde{Q}}[\operatorname{Tr}(\delta^{2})]&\leq\mathbf{E}_{Q}[\operatorname{Tr}(\delta^{2})]=O\left(\frac{LN^{2}}{KM}\right).\end{aligned}

We then estimate

	$\displaystyle d_{\mathrm{TV}}(Q\mathbf{G}Q^{T},\mathbf{G})$	$\displaystyle\leq d_{\mathrm{TV}}(Q\mathbf{G}Q^{T},\tilde{Q}\mathbf{G}\tilde{Q}^{T})+d_{\mathrm{TV}}(\tilde{Q}\mathbf{G}\tilde{Q}^{T},\mathbf{G})$
(2.17)			$\displaystyle\leq\mathbf{P}[\\|\delta(Q)\\|>1/2]+\sqrt{\frac{1}{2}D_{\mathrm{KL}}(\tilde{Q}\mathbf{G}\tilde{Q}^{T}\|\|\mathbf{G})}.$

For any $\tilde{Q}$ which is independent of $\mathbf{G}$ , the density of $Y=\tilde{Q}\mathbf{G}\tilde{Q}^{T}$ can be directly seen to be $\mathbf{E}_{\tilde{Q}}[f_{\tilde{Q}}(Z)]$ . Thus using Jensen’s inequality with $x\mapsto x\log x$ convex, and that $D_{\mathrm{KL}}(\tilde{q}\mathbf{G}\tilde{q}^{T}||\mathbf{G})$ is bounded for example via (2.10),

	$\displaystyle D_{\mathrm{KL}}(\tilde{Q}\mathbf{G}\tilde{Q}^{T}\|\|\mathbf{G})$	$\displaystyle=\int\frac{\mathbf{E}_{\tilde{Q}}[f_{\tilde{Q}}(Z)]}{g(Z)}\log\left[\frac{\mathbf{E}_{\tilde{Q}}[f_{\tilde{Q}}(Z)]}{g(Z)}\right]g(Z)\,dZ$
		$\displaystyle\leq\int\mathbf{E}_{\tilde{Q}}\left[\frac{f_{\tilde{Q}}(Z)}{g(Z)}\log\left[\frac{f_{\tilde{Q}}(Z)}{g(Z)}\right]\right]g(Z)\,dZ$
(2.18)			$\displaystyle=\mathbf{E}_{q\sim{\tilde{Q}}}[D_{\mathrm{KL}}(q\mathbf{G}q^{T}\|\|\mathbf{G})].$

Using that $\tilde{Q}$ is invertible and $\|\delta(\tilde{Q})\|\leq 1/2$ by construction, (2.10) and (2.16) then imply

(2.19)

\displaystyle\mathbf{E}_{\tilde{q}\sim\tilde{Q}}[D_{\mathrm{KL}}(\tilde{q}\mathbf{G}\tilde{q}^{T}||\mathbf{G})]\leq O\left(\frac{LN^{3}}{KM}\right).

Applying this and (2.15) in (2.17) gives

(2.20)

\displaystyle d_{\mathrm{TV}}(Q\mathbf{G}Q^{T},\mathbf{G})

\displaystyle\leq O\left(\sqrt{\frac{(M-K)N^{3}}{KM}}\right),

which with (2.6) implies (2.1). ∎

The proof of Theorem 1.1 for $N=o(\sqrt{K})$ diverges from the above proof by proving a sharper bound on $d_{\mathrm{TV}}(Q\mathbf{G}Q^{T},\mathbf{G})$ , using $\chi^{2}$ -divergence and a Bakry–Émery concentration result [BE85, KM16]. This proof is given in Appendix A. For the rest of the main text of this paper, we will assume the full Theorem 1.1 is proved.

3. Proof of Theorem 1.2

In this section, we prove Theorem 1.2 on closeness of (properly scaled) $GG^{T}$ and $\mathbf{G}\sim\mathcal{G}_{N}^{\mathrm{sym}}$ . It suffices to prove this for $N=O(\sqrt{K})$ . The sparse result in [SMG26, Theorem 1.5] shows that for $N\leq K$ and $NK=o(M)$ , that

(3.1)

\displaystyle d_{\mathrm{TV}}(MU_{NK}U_{NK}^{T},G_{NK}G_{NK}^{T})

\displaystyle\leq O\left(\sqrt{\frac{NK}{M}}\right),

where $G_{NK}$ is an $N\times K$ matrix of iid complex standard Gaussians. We use this result with Theorem 1.1 (with the full $N=O(\sqrt{K})$ result) to prove Theorem 1.2. Choose e.g. $M=K^{3}$ (any $\omega(K^{2})$ will do). Then for $N=O(K^{1/2})$ ,

	$\displaystyle d_{\mathrm{TV}}(\mathbf{G},G_{NK}G_{NK}^{T}/\sqrt{K})$	$\displaystyle\leq d_{\mathrm{TV}}(\mathbf{G},MK^{-1/2}U_{NK}U_{NK}^{T})+d_{\mathrm{TV}}(MU_{NK}U_{NK}^{T},G_{NK}G_{NK}^{T})$
		$\displaystyle\leq O\left(\frac{N}{\sqrt{K}}\right)+O\left(\sqrt{\frac{NK}{M}}\right)\leq O\left(\frac{N}{\sqrt{K}}\right),$

as desired. ∎

4. Hardness argument

In this section, we prove Theorem 1.4. There are two main differences from the boson sampling argument of [AA13]:

(1)

We have a choice of squeezing parameter $s$ and a variable photon number $N$ .
(2)

We lack a general instance generating/rejection sampling lemma to hide $X\in\mathcal{G}_{N}^{\mathrm{sym}}$ exactly as a submatrix of $UI_{K}U^{T}$ .

For (1), because the output probabilities (1.1) involve the squeezing parameter $s$ , the additive error threshold obtained will depend on the choice of $s$ . We will choose $s$ to minimize this additive error, which will allow us to state the $|\mathrm{GHE}|_{\pm}^{2}$ problem in terms of additive error $\varepsilon\cdot\mathbf{E}|\operatorname{Haf}(\mathbf{G})|^{2}$ . Difference (2) will be resolved by using an approximate instead of exact instance generating method in FPostBPP. Since $\text{PostBPP}\subseteq\text{BPP}^{\text{NP}}$ [AA13, §2], this will not change the end complexity class. Proposition 1.3 shows why one cannot obtain an exact instance generating lemma via rejection sampling as in [AA13, Lemma 5.7] for Gaussian boson sampling with $K=\alpha M$ , $0<\alpha<1/2$ . We note that previous GBS hardness arguments avoided the issue in (2) by either working in the sparse squeezer regime [HKS+17, KHS+19], where instance generating follows from instance generating for Fock boson sampling [AA13], or by considering $K=M$ [SMG26] in which case the bound to use the rejection sampling of [AA13] is easily found to hold, or by considering the variant bipartite Gaussian boson sampling [GBA+22, BBD+26], which uses a different set-up with two-mode squeezed states with parameters tuned to implement a specific matrix. For the setting here with arbitrary $K$ including $K=\alpha M$ , $0<\alpha<1/2$ , we will implement an approximate instance generating procedure.

For difference (1), note that for an initial state consisting of $K$ single-mode squeezed states with equal squeezing parameters $s_{i}=s>0$ , and any photon count $N\in 2\mathbb{N}$ [KHS+19],

(4.1)

\displaystyle\mathbf{P}[N]

\displaystyle=\binom{N/2+K/2-1}{N/2}\frac{\tanh^{N}s}{\cosh^{K}s},

where the binomial coefficient is a generalized binomial coefficient allowing for non-integer $K/2$ . The average number of photons for such an input state is $K\sinh^{2}s$ . We will have $N=O(\sqrt{K})$ , and will choose squeezing $s$ (to $\operatorname{poly}(N,1/\delta,1/\varepsilon)$ bits) so that, up to exponentially small error,

(4.2)

\displaystyle N=K\sinh^{2}s,\quad\text{which implies}\quad\frac{\cosh^{K}s}{\tanh^{N}s}=\frac{\left(1+\frac{N}{K}\right)^{K/2}}{\left(\frac{N}{N+K}\right)^{N/2}}=\frac{e^{N/2}K^{N/2}}{N^{N/2}}\exp\left[\frac{N^{2}}{4K}+o(1)\right],

using that $N^{2}=O(K)$ . This choice of $s$ can be seen to minimize the ratio $\frac{\cosh^{K}s}{\tanh^{N}s}$ (e.g. by taking logarithmic derivative), which will appear in the additive error threshold (4.25). Also, note that for $N=o(\sqrt{K})$ , as $N\to\infty$ ,

(4.3)

\displaystyle\mathbf{P}[N]

\displaystyle=\frac{\Gamma(N/2+K/2)}{(N/2)!\Gamma(K/2)}\frac{\left(\frac{N}{N+K}\right)^{N/2}}{\left(1+\frac{N}{K}\right)^{K/2}}=\frac{1+o(1)}{\sqrt{\pi N}}.

Since this is no worse than inverse polynomial, we could postselect on observing exactly $N$ photons, at the cost of repeating the GBS experiment polynomial-many more times. However, this postselection will actually be unnecessary in the hardness argument, since we will be using Stockmeyer’s algorithm [Sto85] to estimate output probabilities. The main point of choosing $s$ as in (4.2) is to minimize the additive error threshold in (4.25). We do not use postselection on $N$ or high collision-free probability anywhere in the hardness reduction.

In order to give a hardness reduction, we need to consider finite-precision inputs and outputs to a Turing machine. In particular, the Gaussian and unitary matrices involve real and complex numbers, which must be rounded to finite-precision. As in [AA13, Footnote 19], we encode each entry of a unitary matrix in binary to $\operatorname{poly}(M,1/\varepsilon,1/\delta)$ bits. In general, the resulting binary description $v$ is not actually unitary, but we can associate each description $v$ to a unitary $U(v)\in\mathrm{U}(M)$ obtained by Gram–Schmidt orthonormalization of the columns with typically very small error [Aar03, Lemma 7.2]. We let $\mathrm{U}^{\prime}(M)$ denote the set of such finite binary descriptions $v$ , and let $U(v)$ denote its associated unitary $U\in\mathrm{U}(M)$ . Based on the above, we make the following finite-precision implementation assumption, which we expect can be proved via careful accounting with Gram–Schmidt or QR factorization. We note that Assumption 1(a) is one possible route, though not the only one, to rigorously implement the finite-precision vs Haar sampling also used in [AA13].

Assumption 1.

Let $\eta>0$ be a finite precision and $\rho>0$ a failure probability. There is a bit precision $Q=\operatorname{poly}(M,\log(1/\eta),\log(1/\rho))$ , a number of random bits $R=\operatorname{poly}(M,Q)$ , and a polynomial-time algorithm which maps $\{0,1\}^{R}\ni\xi\mapsto v=v_{\xi}$ for $v$ a finite binary description with associated $M\times M$ unitary $U(v)$ as described above, such that the following hold:

(1)

There is a coupling $(v,H)$ , with $H\sim\mathcal{H}_{M}$ Haar random, such that for $\xi$ uniform,

$\displaystyle\mathbf{P}[\|U(v)-H\|_{\mathrm{op}}>\eta/4]\leq\rho.$
(2)

From $(v,S)$ , where $S$ is any subset $S\subset\{1,\ldots,M\}$ , one can compute a finite-precision binary matrix $\hat{Z}(v,S)$ in deterministic polynomial time such that

$\displaystyle\|\hat{Z}(v,S)-(U(v)I_{K}U(v)^{T})_{S,S}\|_{\infty}\leq\eta/2.$

As a result, combining (a) and (b), for this coupling and for $S$ an independent uniform random size $N$ subset,

(4.4)

\displaystyle\mathbf{P}[\|\hat{Z}(v,S)-(HI_{K}H^{T})_{S,S}\|_{\infty}>\eta]\leq\rho.

We let $\mu_{Q,M}$ denote the law of $v_{\xi}$ for uniform $\xi$ .

The main property we need from Assumption 1 is (4.4), which says we can produce a binary finite-precision matrix $\hat{Z}(v,S)$ which is generally a good entry-wise approximation to $(HI_{K}H^{T})_{S,S}$ . We need the finite-precision matrix so we can apply an NP-oracle in the proof of Lemma 4.1 below. The failure probability $\rho$ in (4.4) accounts for ill-conditioned starting Gaussian matrices which may have poor finite-precision approximation during the coupling procedure.

Compared to [AA13], which generally referred directly to the continuum unitary and Gaussian distributions after the preliminary finite-precision discussions, we will need to be more precise, and will keep track of the finite-precision rounding. We will need to use the complexity class PostBPP (also called $\text{BPP}_{\text{path}}$ ), which is BPP with postselection, defined precisely in e.g. [AA13, §2]. The main property we need for this class is $\text{PostBPP}\subseteq\text{BPP}^{\text{NP}}$ [HHT97, BGP00]. We will also use the function/search analogue FPostBPP. To resolve the difference (2) from above, we prove (under Assumption 1)

Lemma 4.1 (approximate instance generation).

Consider $X\sim\mathcal{G}_{N}^{\mathrm{sym}}$ , and let $\delta,\varepsilon$ be error parameters. Let $K=\Omega(N^{2}/\delta^{2})$ , $M=O(\operatorname{poly}(N,1/\delta))$ , $N\leq K\leq M$ , and $[M]:=\{1,\ldots,M\}$ . Choose a scale $\gamma=2^{-p}$ for some $p=\operatorname{poly}(N,1/\delta,1/\varepsilon)$ . Then there is an $\text{FBPP}^{\text{NP}}$ algorithm $\mathcal{A}$ (running in $\operatorname{poly}(N,1/\delta,1/\varepsilon)$ time) which, on a sufficiently fine $\operatorname{poly}(N,1/\delta,1/\varepsilon)$ -bit precision implementation, takes as input a matrix $X\sim\mathcal{G}_{N}^{\mathrm{sym}}$ , and outputs either $\bot$ (failure) or a pair $(v,S^{*})\in\mathrm{U}^{\prime}(M)\times\{\text{size $N$ subsets of $[M]$}\}$ . It succeeds with probability $\geq 1-O(\delta)$ over $X$ and $\mathcal{A}$ . Conditioned on succeeding, the output $(v,S^{*})$ satisfies

(1)

$\|(U(v)I_{K}U(v)^{T})_{S^{*},S^{*}}-M^{-1}K^{1/2}X\|_{\infty}\leq O(\gamma)$ , where $\|\cdot\|_{\infty}$ denotes the $\ell^{\infty}$ maximum entrywise norm, and $U(v)$ the unitary associated with the finite-precision description $v$ .
(2)

Let $\mathcal{L}(v,S^{*})$ be the law of $(v,S^{*})$ averaged over $X$ in successful trials, and let $\mu_{Q,M}\otimes\operatorname{Unif}$ denote the product measure of $\mu_{Q,M}$ on $\mathrm{U}^{\prime}(M)$ , and the uniform distribution on size $N$ subsets of $[M]$ . Then

(4.5) $\displaystyle d_{\mathrm{TV}}(\mathcal{L}(v,S^{*}),\mu_{Q,M}\otimes\operatorname{Unif})\leq O(\delta).$

Proof of Lemma 4.1.

Consider the true distributions $X^{\prime}=M^{-1}K^{1/2}X$ and $H\sim\mathcal{H}_{M}$ , and let $Z:=(HI_{K}H^{T})_{S,S}$ for a uniformly random size $N$ subset $S\subset[M]$ . Let $R_{\gamma}$ denote rounding the real and imaginary part of every entry of a matrix to the interval $[j\gamma,(j+1)\gamma)$ , $j\in\mathbb{Z}$ , containing it. The overall idea, ignoring finite-precision implementation for now, is: Given $X^{\prime}$ , we want to generate random $(H,S)$ and postselect on the event $R_{\gamma}(Z)=R_{\gamma}(X^{\prime})$ , i.e. up to $\gamma$ error, $X^{\prime}$ appears as the submatrix $Z=(HI_{K}H^{T})_{S,S}$ . Intuitively, this is a postselection problem so can be done in $\text{FPostBPP}\subseteq\text{FBPP}^{\text{NP}}$ . However, we need to consider the finite-precision implementation details to properly apply the NP oracle in $\text{FBPP}^{\text{NP}}$ [BGP00]. We do this as follows.

(1)

Truncation. First note we can restrict to $X$ with maximum entry size $\|X\|_{\infty}\leq C_{1}\sqrt{\log(N/\delta)}$ for large enough $C_{1}$ , as this occurs with high probability $1-O(\delta)$ by Gaussian suprema and concentration bounds.⁶⁶ 6 If $X_{t}$ is $\sigma^{2}$ -subgaussian for every $t\in T$ , then $\mathbf{P}[\sup_{t\in T}X_{t}\geq\sqrt{2\sigma^{2}\log|T|}+x]\leq e^{-x^{2}/2\sigma^{2}}$ ; see e.g. [vH16, §5]. In this case, $|T|=N(N+1)/2$ . Then let $X^{\prime}=M^{-1}K^{1/2}X$ conditioned on $\|X\|_{\infty}\leq C_{1}\sqrt{\log(N/\delta)}$ . The cutoff prevents arbitrarily large entries, and allows for the required finite-precision implementation. Note that $d_{\mathrm{TV}}(Z,X^{\prime})\leq d_{\mathrm{TV}}(Z,M^{-1}K^{1/2}X)+d_{\mathrm{TV}}(M^{-1}K^{1/2}X,X^{\prime})\leq O(\delta)$ , using Theorem 1.1, invariance of Haar measure under row permutations, and the TVD coupling inequality. (If $X^{\prime}\overset{d}{=}Y|\mathcal{E}$ , then by a coupling argument, $d_{\mathrm{TV}}(X^{\prime},Y)\leq\mathbf{P}[\mathcal{E}^{c}]$ .)

(2)

Finite-precision cells and bounds. Choose finite precision $\eta=2^{-q}$ with say $\eta\leq\min(\delta\gamma^{2},\delta 2^{-N})$ , which requires only $q\sim\operatorname{poly}(N,1/\delta,1/\varepsilon)$ bits of precision. Taking $\rho=O(\delta)$ in (4.4) of Assumption 1, with probability $1-O(\delta)$ we can produce a finite-precision matrix $\hat{Z}(v,S)$ such that $\|\hat{Z}(v,S)-Z\|_{\infty}\leq\eta$ . There is also a finite-precision matrix $\hat{X}^{\prime}$ such that $\|\hat{X}^{\prime}-X^{\prime}\|_{\infty}\leq\eta$ . We want to show that just like $d_{\mathrm{TV}}(R_{\gamma}(Z),R_{\gamma}(X^{\prime}))\leq d_{\mathrm{TV}}(Z,X^{\prime})=O(\delta)$ for the continuum distributions, that

(4.6)

\displaystyle d_{\mathrm{TV}}(R_{\gamma}(\hat{Z}(v,S)),R_{\gamma}(\hat{X}^{\prime}))=O(\delta).

This bound will ensure the language to sample from later in (4.10) is non-empty with high probability, and that the resulting distribution of $(v,S^{*})$ satisfies part (ii) of the lemma.

We start by applying the triangle inequality, giving

(4.7)

d_{\mathrm{TV}}(R_{\gamma}(\hat{Z}(v,S)),R_{\gamma}(\hat{X}^{\prime}))\\ \leq d_{\mathrm{TV}}(R_{\gamma}(\hat{Z}(v,S)),R_{\gamma}(Z))+d_{\mathrm{TV}}(R_{\gamma}(Z),R_{\gamma}(X^{\prime}))+d_{\mathrm{TV}}(R_{\gamma}(X^{\prime}),R_{\gamma}(\hat{X}^{\prime})).

To bound the first and third distances on the right side, we need to bound the probability that $\hat{Z}(v,S)$ or $\hat{X}^{\prime}$ may have moved $\gamma$ -cells during the finite-precision implementation/rounding process. Intuitively, because we can take much finer precision $\eta$ than $\gamma$ , this probability is very small. More formally, let $\partial_{\eta}=\partial_{\eta,\gamma}$ denote the set of matrices with some real or imaginary coordinate within $\eta$ of the grid $\gamma\mathbb{Z}$ . Note that for a matrix outside of $\partial_{\eta}$ , moving by $\eta$ cannot move across the coarser $\gamma$ -cell boundaries. The boundary estimate Lemma 4.2 below then gives $\mathbf{P}[X^{\prime}\in\partial_{\eta}]=O(\delta)$ , and

\displaystyle\mathbf{P}[Z\in\partial_{\eta}]

\displaystyle\leq\mathbf{P}[X^{\prime}\in\partial_{\eta}]+d_{\mathrm{TV}}(Z,X^{\prime})=O(\delta).

If the coupling in Assumption 1 is succesful, then $\|\hat{Z}(v,S)-Z\|_{\infty}\leq\eta$ . Thus the rounded cells $R_{\gamma}(\hat{Z}(v,S))$ and $R_{\gamma}(Z)$ can differ only if the coupling in Assumption 1 failed, or if $Z\in\partial_{\eta}$ . Thus by the TVD coupling bound and Assumption 1,

(4.8)

\displaystyle d_{\mathrm{TV}}(R_{\gamma}(\hat{Z}(v,S)),R_{\gamma}(Z))

\displaystyle\leq\rho+\mathbf{P}[Z\in\partial_{\eta}]=O(\delta).

Also by TVD coupling lemma,

(4.9)

\displaystyle d_{\mathrm{TV}}(R_{\gamma}(X^{\prime}),R_{\gamma}(\hat{X}^{\prime}))

\displaystyle\leq\mathbf{P}[X^{\prime}\in\partial_{\eta}]=O(\delta).

The remaining term in the middle of the right side of (4.7) is bounded using the hiding property Theorem 1.1 and TVD conditioning property, as $d_{\mathrm{TV}}(R_{\gamma}(Z),R_{\gamma}(X^{\prime}))\leq d_{\mathrm{TV}}(Z,X^{\prime})=O(\delta)$ . This gives (4.6).

(3)

Uniform NP witness generation with an NP oracle. The precise finite-precision problem is now: Given the cell $R_{\gamma}(\hat{X}^{\prime})$ , post-select on the event $R_{\gamma}(\hat{Z}(v,S))=R_{\gamma}(\hat{X}^{\prime})$ . Recall we view $\hat{Z}(v,S)=\hat{Z}(v_{r},S_{t})$ as generated via a random-string model $r,t\in\{0,1\}^{\operatorname{poly}(N,1/\delta,1/\varepsilon)}\mapsto(v_{r},S_{t})$ . Note, in order to sample uniformly from the $\binom{M}{N}$ choices of $S$ , which need not divide $2^{L}$ , we actually consider uniform random $t\in\{0,1\}^{\ell}$ for $\ell=\left\lceil\log_{2}\binom{M}{N}\right\rceil$ , and accept if $t<\binom{M}{N}$ , then map each of those $t$ to a subset $S=S_{t}$ . For the above finite-precision problem, we want to sample a uniform $(r,t)$ from the language

(4.10)

\displaystyle L_{\hat{X}^{\prime}}=\{r\in\{0,1\}^{\operatorname{poly}(N,1/\delta,1/\varepsilon)},t\in\{0,1\}^{\ell}:R_{\gamma}(\hat{Z}(v_{r},S_{t}))=R_{\gamma}(\hat{X}^{\prime})\text{ and }t<\binom{M}{N}\},

and then output $(v_{r},S_{t})$ . Due to the TVD bound (4.6), and considering the set of $\gamma$ -cells $\{c:\text{not attainable by }R_{\gamma}(\hat{Z}(v_{r},S_{t})))\}$ , we see $L_{\hat{X}^{\prime}}$ is nonempty with probability $1-O(\delta)$ over $\hat{X}^{\prime}$ . When $L_{\hat{X}^{\prime}}$ is nonempty, sampling a uniform pair $(r,t)$ can be done in probabilistic polynomial time with an NP oracle via [BGP00]: the language $L_{\hat{X}^{\prime}}$ is in NP (and P), so if $L_{\hat{X}^{\prime}}$ is nonempty, then [BGP00] generates a uniform random witness $(r,t)\in L_{\hat{X}^{\prime}}$ with success probability at least $0.2$ . Additionally, we can obtain success probability $1-O(\delta)$ using standard amplification with $O(\log\delta^{-1})$ trials in the FPostBPP to $\text{FBPP}^{\text{NP}}$ implementation. In total, conditioned on successful output (including $\mathcal{L}_{\hat{X}^{\prime}}$ nonempty, which can be checked by seeing failure or verifying if the output is in $\mathcal{L}_{\hat{X}^{\prime}}$ ), we obtain a uniform random sample $(v,S^{*})$ from the distribution $\mu_{Q,M}\otimes\operatorname{Unif}$ conditioned on the event $R_{\gamma}(\hat{Z}(v,S^{*}))=R_{\gamma}(\hat{X}^{\prime})$ , which is generated in $\text{FBPP}^{\text{NP}}$ in time $\operatorname{poly}(N,1/\delta,1/\varepsilon)$ and with success probability $1-O(\delta)$ .

Since $\|\hat{Z}(v,S^{*})-(U(v)I_{K}U(v)^{T})_{S^{*},S^{*}}\|_{\infty}\leq\eta/2$ and $\|X^{\prime}-\hat{X}^{\prime}\|_{\infty}\leq\eta$ , we obtain

(4.11)

\displaystyle\|(U(v)I_{K}U(v)^{T})_{S^{*},S^{*}}-X^{\prime}\|_{\infty}

\displaystyle\leq O(\gamma)+O(\eta)=O(\gamma).

Since we can include $\|X\|_{\infty}>C_{1}\sqrt{\log(N/\delta)}$ in the failure probability, (4.11) gives part (i) of the lemma.

(4)

It remains to prove (ii) of the lemma. Given a cell $c=R_{\gamma}(\hat{X}^{\prime})$ attainable by $R_{\gamma}(\hat{Z}(v_{r},S_{t}))$ , successful outputs $(v,S^{*})$ are distributed as $\mu_{Q,M}\otimes\operatorname{Unif}$ conditioned on $R_{\gamma}(\hat{Z}(w,S))=c$ ; call this conditional measure $K_{c}$ , for $c$ attainable by $R_{\gamma}(\hat{Z}(v,S))$ . The total distribution of $(v,S^{*})$ is then $\sum_{c}\mathbf{P}[R_{\gamma}(\hat{X}^{\prime})=c|\text{success}]K_{c}$ . We can also decompose $\mu_{Q,M}\otimes\operatorname{Unif}$ as $\sum_{c}\mathbf{P}[R_{\gamma}(\hat{Z}(w,S))=c]K_{c}$ for $(w,S)\sim\mu_{Q,M}\otimes\operatorname{Unif}$ .

Using the characterizations (1.3) and (1.4) of TVD, we thus obtain

	$\displaystyle d_{\mathrm{TV}}(\mathcal{L}(v,S^{*}),\mu_{Q,M}\otimes\operatorname{Unif})$	$\displaystyle=\sup_{f:\\|f\\|_{\infty}\leq 1}\frac{1}{2}\sum_{c}(\mathbf{P}[R_{\gamma}(\hat{X}^{\prime})=c\|\text{success}]-\mathbf{P}[R_{\gamma}(\hat{Z}(w,S))=c])K_{c}(f)$
		$\displaystyle\leq d_{\mathrm{TV}}(R_{\gamma}(\hat{X}^{\prime})\|\text{success},R_{\gamma}(\hat{Z}(w,S)))$
		$\displaystyle\leq d_{\mathrm{TV}}(R_{\gamma}(\hat{X}^{\prime}),R_{\gamma}(\hat{Z}(w,S)))+O(\delta)=O(\delta),$

using (4.6) and that $\mathbf{P}[\text{failure}]=O(\delta)$ . This completes the proof of Lemma 4.1(ii).

∎

The following lemma was used in the proof of Lemma 4.1.

Lemma 4.2 (boundary estimate).

Let $\eta\leq\gamma/2$ , and let $\partial_{\eta}=\partial_{\eta,\gamma}$ denote the set of matrices with some real or imaginary coordinate within $\eta$ of the grid $\gamma\mathbb{Z}$ . Then for $X^{\prime}=M^{-1}K^{1/2}X$ with $X\sim\mathcal{G}_{N}^{\mathrm{sym}}$ ,

(4.12)

\displaystyle\mathbf{P}[X^{\prime}\in\partial_{\eta}]

\displaystyle\leq CN^{2}\eta(\gamma^{-1}+MK^{-1/2}).

If we consider $X^{\prime}$ conditioned on a probability $1-\epsilon$ event, then the above bound holds with an additional $+\epsilon$ term.

Proof.

We do a union bound over the $O(N^{2})$ independent entries of $X$ . For a single real Gaussian $Y\sim\mathcal{N}(0,\sigma^{2})$ with density function $f_{Y}$ , which is decreasing for $y\geq 0$ ,

	$\displaystyle\mathbf{P}[\operatorname{dist}(Y,\gamma\mathbb{Z})\leq\eta]$	$\displaystyle=\sum_{j\in\mathbb{Z}}\int_{j\gamma-\eta}^{j\gamma+\eta}f_{Y}(y)\,dy$
		$\displaystyle\leq 2\eta\\|f_{Y}\\|_{\infty}+2\sum_{j=1}^{\infty}2\eta f_{Y}(j\gamma-\eta)$
(4.13)			$\displaystyle\leq\frac{2\eta}{\sigma\sqrt{2\pi}}+\frac{4\eta}{\gamma-\eta}\sum_{j=1}^{\infty}\int_{j\gamma-\gamma}^{j\gamma-\eta}f_{Y}(y)\,dy\leq C\eta(\sigma^{-1}+\gamma^{-1}).$

Since entries of $X^{\prime}$ have standard deviation $\sigma=\Theta(M^{-1}K^{1/2})$ , a union bound gives (4.12). The conditioning statement holds since for a random variable $Z$ and event $\mathcal{E}$ , letting $Y:\overset{d}{=}Z|\mathcal{E}$ , then $d_{\mathrm{TV}}(Z,Y)\leq\mathbf{P}[\mathcal{E}^{c}]$ , by using the coupling bound. (Set $Y=Z$ on $\mathcal{E}$ , and $Y\sim\mathcal{L}(Z|\mathcal{E})$ on $\mathcal{E}^{c}$ .) ∎

Because we round to $\operatorname{poly}(N,1/\delta,1/\varepsilon)$ bits, we want to check the resulting hafnian is not changed too much. Suppose $\|X\|_{\infty},\|Y\|_{\infty}\leq B$ . Then for $N\in 2\mathbb{N}$ ,

	$\displaystyle\left\|\|\operatorname{Haf}X\|^{2}-\|\operatorname{Haf}Y\|^{2}\right\|$	$\displaystyle=\Bigg\|\sum_{\pi,\pi^{\prime}\in\mathcal{P}_{2}(N)}\left(\prod_{\{i,j\}\in\pi}X_{ij}\prod_{\{i^{\prime},j^{\prime}\}\in\pi^{\prime}}\bar{X}_{i^{\prime}j^{\prime}}-\prod_{\{i,j\}\in\pi}Y_{ij}\prod_{\{i^{\prime},j^{\prime}\}\in\pi^{\prime}}\bar{Y}_{i^{\prime}j^{\prime}}\right)\Bigg\|$
		$\displaystyle\leq\sum_{\pi,\pi^{\prime}\in\mathcal{P}_{2}(N)}NB^{N-1}\\|X-Y\\|_{\infty}$
(4.14)			$\displaystyle\leq(N-1)!!^{2}NB^{N-1}\\|X-Y\\|_{\infty}.$

The rounding process changes a matrix by at most $\|X-Y\|_{\infty}\leq\gamma=2^{-\operatorname{poly}(N,1/\delta,1/\varepsilon)}$ . Since $B=O(M^{-1}K^{1/2}\sqrt{\log(N/\delta)})$ for the sampled $U$ in Lemma 4.1, $\gamma$ is easily chosen so that (4.14) is $\leq 2^{-\operatorname{poly}(N,1/\delta,1/\varepsilon)}$ for such matrices.

Proof of Theorem 1.4.

Suppose we had an oracle $\mathcal{O}$ for limited instances of approximate Gaussian boson sampling as described in the hypotheses. $\mathcal{O}$ takes as input a string $r\in\{0,1\}^{\operatorname{poly}(M,1/\delta,1/\varepsilon)}$ , a finite binary description $v$ representing an $M\times M$ unitary matrix $U=U(v)$ , a squeezing parameter $s>0$ for the first $K=K(M)$ input modes, and an error bound $\varepsilon>0$ . We may write the Gaussian boson sampler description as $A=A(v,K,s)$ . Over uniform random $r$ (representing randomness in the Gaussian boson sampling output), $\mathcal{O}$ outputs a distribution $\mathcal{D}_{\mathcal{O}}(A,\varepsilon)$ which is $\varepsilon$ -close to $\mathcal{D}_{A}$ , the exact Gaussian boson sampling output distribution for linear optical unitary $U(v)$ .

Let $X\sim\mathcal{G}_{N}^{\mathrm{sym}}$ be an input matrix and let $\varepsilon,\delta>0$ be the given error parameters. Note that it suffices to solve $|\mathrm{GHE}|_{\pm}^{2}$ with probability at least $1-O(\delta)$ , since we can simply rescale $\delta$ to get the probability to $1-\delta$ . So we want to approximate $|\operatorname{Haf}(X)|^{2}$ to within additive error $\varepsilon\cdot(N-1)!!$ with success probability at least $1-O(\delta)$ over $X$ .

Choose $(K(M),M)$ with $M=\operatorname{poly}(K(M))$ , and squeezing parameter $s$ , so that

(4.15)

\displaystyle N\leq c_{1}\delta K^{1/2},\qquad N=K\sinh^{2}s,

for a constant $c_{1}$ chosen so that the TVD error bound in (4.5) is e.g. $\leq\delta/8$ . Rescaling $X^{\prime}:=M^{-1}K^{1/2}X$ , we then want to approximate $|\operatorname{Haf}(X^{\prime})|^{2}$ to within additive error $\varepsilon\cdot(N-1)!!(M^{-1}K^{1/2})^{N}$ .

By Lemma 4.1, for a sufficiently fine numerical precision, with probability $1-O(\delta)$ , we can efficiently generate a random finite binary description $v=v_{\xi}$ of an $M\times M$ unitary $U=U(v)$ , and a size $N$ subset $S^{*}$ of $[M]$ , such that

(4.16)

\displaystyle\|(UI_{K}U^{T})_{S^{*},S^{*}}-X^{\prime}\|_{\infty}\leq 2^{-\operatorname{poly}(N,1/\delta,1/\varepsilon)},\quad\text{and}\quad d_{\mathrm{TV}}(\mathcal{L}(v,S^{*}),\mu_{Q,M}\otimes\operatorname{Unif})\leq O(\delta).

Recall $\mu_{Q,M}$ is the measure on the finite binary descriptions $\mathrm{U}^{\prime}(M)$ defined in Assumption 1.

Conditioned on success of the above instance generating procedure, take $v$ and send $A=A(v,K,s)$ to the GBS oracle $\mathcal{O}$ . We will also allow an adversary to know $N$ , since it will not affect the argument, and since they could possibly already guess a narrow range of $N$ we are interested in based on the squeezing parameters provided. Let $\beta$ be an error bound which is polynomial in $\varepsilon$ and $\delta$ . Then more formally, denoting $0^{1/\beta}=0\cdots 0$ with $\lceil 1/\beta\rceil$ zeros, we send the input $\langle A(v,K,s),0^{1/\beta},r\rangle$ , for $r\in\{0,1\}^{\operatorname{poly}(M,1/\delta,1/\varepsilon)}$ a random string, to the oracle $\mathcal{O}$ . As $r$ is varied, $\mathcal{O}$ returns a sample from $\mathcal{D}_{A}^{\prime}$ with $\|\mathcal{D}_{A}-\mathcal{D}_{A}^{\prime}\|_{\mathrm{TV}}\leq\beta$ .

For collision-free $S$ with $N$ photon counts, let

\displaystyle p_{S}(v):=\mathbf{P}_{\mathcal{D}_{A}}[S],\quad q_{S}(v):=\mathbf{P}_{\mathcal{D}_{A}^{\prime}}[S]=\mathbf{P}_{r}[\mathcal{O}(A(v,K,s),0^{1/\beta},r)=S].

The probability $q_{S^{*}}$ will be approximated in $\text{FBPP}^{\text{NP}^{\mathcal{O}}}$ by Stockmeyer’s algorithm [Sto85] as usual. The probability $p_{S^{*}}$ is

(4.17)

\displaystyle p_{S^{*}}

\displaystyle=\frac{\tanh^{N}(s)}{\cosh^{K}(s)}|\operatorname{Haf}[(UI_{K}U^{T})_{S^{*},S^{*}}]|^{2}.

Given $p_{S^{*}}$ or a good approximation to $p_{S^{*}}$ , this with (4.14) and Lemma 4.1(i) will then let us estimate $|\operatorname{Haf}(X^{\prime})|^{2}$ .

We want to show that $p_{S^{*}}$ and $q_{S^{*}}$ are close with high probability over $X$ and $v$ . We can average over the $\binom{M}{N}$ possible $N$ -photon collision-free outputs $S$ as in [AA13, §5.2] to obtain

(4.18)

\displaystyle\mathbf{E}_{|S|=N}[|p_{S}-q_{S}|]

\displaystyle\leq\frac{\sum_{|S|=N}|p_{S}-q_{S}|}{\binom{M}{N}}\leq\frac{2\|\mathcal{D}_{A}-\mathcal{D}^{\prime}_{A}\|_{\mathrm{TV}}}{\binom{M}{N}}<C_{1}\beta\frac{N!}{M^{N}},

using that $M\geq K\geq c_{1}^{-2}\delta^{-2}N^{2}$ . Then taking $\beta=\varepsilon\delta$ , Markov’s inequality gives

(4.19)

\displaystyle\mathbf{P}_{|S|=N}\left[|p_{S}-q_{S}|>\frac{\varepsilon}{2}\frac{N!}{M^{N}}\right]

\displaystyle\leq\frac{O(\beta)}{\varepsilon}\leq O(\delta).

We will use Lemma 4.1(ii) to obtain a bound on $|p_{S^{*}}-q_{S^{*}}|$ . First, if $(w,S)\sim\mu_{Q,M}\otimes\operatorname{Unif}$ , then since $S$ is uniform, the above implies

(4.20)

\displaystyle\mathbf{P}_{(w,S)}\left[|p_{S}(w)-q_{S}(w)|>\frac{\varepsilon}{2}\frac{N!}{M^{N}}\right]

\displaystyle\leq O(\delta).

Applying Lemma 4.1(ii), we have

	$\displaystyle\mathbf{P}_{(v,S^{})}\left[\|p_{S^{}}(v)-q_{S^{*}}(v)\|>\frac{\varepsilon}{2}\frac{N!}{M^{N}}\right]$	$\displaystyle\leq\mathbf{P}_{(w,S)}\left[\|p_{S}(w)-q_{S}(w)\|>\frac{\varepsilon}{2}\frac{N!}{M^{N}}\right]+O(\delta)$
(4.21)			$\displaystyle\leq O(\delta).$

Next, as in [AA13, §5.2], we use Stockmeyer’s algorithm to approximate $q_{S^{*}}$ . For any $\alpha>0$ , Stockmeyer’s algorithm [Sto85] applied with $f(r):=\mathbf{1}[\mathcal{O}(A,0^{1/\beta},r)=S^{*}]$ , and standard output probability amplification, gives an estimate $\tilde{q}_{S^{*}}$ in $\operatorname{poly}(M,1/\alpha)$ time such that

(4.22)

\displaystyle\mathbf{P}[|\tilde{q}_{S^{*}}-q_{S^{*}}|>\alpha q_{S^{*}}]<\frac{1}{2^{M}}.

Also, by Markov’s inequality with $\mathbf{E}_{|S|=N}[q_{S}]\leq\binom{M}{N}^{-1}\leq\frac{2N!}{M^{N}}$ , and using the TVD bound in Lemma 4.1(ii) like in (4.21), we have for any $j>1$ ,

(4.23)

\displaystyle\mathbf{P}_{(v,S^{*})}\left[q_{S^{*}}>2j\cdot\frac{N!}{M^{N}}\right]

\displaystyle\leq\frac{1}{j}+O(\delta).

Thus taking $j=4/\delta$ and $\alpha=\varepsilon\delta/16$ , and also using (4.21), we get

	$\displaystyle\mathbf{P}\left[\|\tilde{q}_{S^{}}-p_{S^{}}\|>\varepsilon\frac{N!}{M^{N}}\right]$	$\displaystyle\leq\mathbf{P}\left[\|\tilde{q}_{S^{}}-q_{S^{}}\|>\frac{\varepsilon}{2}\frac{N!}{M^{N}}\right]+\mathbf{P}\left[\|q_{S^{}}-p_{S^{}}\|>\frac{\varepsilon}{2}\frac{N!}{M^{N}}\right]$
		$\displaystyle\leq\mathbf{P}\left[\|\tilde{q}_{S^{}}-q_{S^{}}\|>\alpha q_{S^{}}\right]+\mathbf{P}\left[q_{S^{}}>2j\frac{N!}{M^{N}}\right]+O(\delta)$
(4.24)			$\displaystyle\leq\frac{1}{2^{M}}+\frac{\delta}{4}+O(\delta).$

Since $M\geq K\geq\Omega(N^{2}/\delta^{2})$ , we get $2^{-M}=O(\delta)$ . We combine this with the $O(\delta)$ failure probability from Lemma 4.1. In total, using the approximation $\tilde{q}_{S^{*}}$ , with probability at least $1-O(\delta)$ we can approximate $|\operatorname{Haf}[(UI_{K}U^{T})_{S^{*},S^{*}}]|^{2}$ , and thus also $|\operatorname{Haf}(X^{\prime})|^{2}$ , to additive error

(4.25)

\displaystyle\varepsilon\cdot\frac{\cosh^{K}(s)}{\tanh^{N}(s)}\frac{N!}{M^{N}}+O(2^{-\operatorname{poly}(N,1/\delta,1/\varepsilon)})

\displaystyle\leq\varepsilon\frac{K^{N/2}}{M^{N}}\frac{N^{N/2}}{e^{N/2}}\sqrt{2\pi N}(1+O(\delta)+O(1/N)),

recalling we chose $s$ in (4.2) up to exponentially small error, and have a small rounding error from (4.14). Comparing to

\displaystyle\varepsilon\cdot\mathbf{E}|\operatorname{Haf}(X^{\prime})|^{2}=\varepsilon\frac{(N-1)!!K^{N/2}}{M^{N}}

\displaystyle=\varepsilon\frac{K^{N/2}}{M^{N}}\frac{N^{N/2}}{e^{N/2}}\sqrt{2}(1+o(1)),

we see (4.25) differs by a factor of $\sqrt{\pi N}$ . But $\varepsilon$ can absorb $\operatorname{poly}(N)$ factors (just run the procedure with $\varepsilon^{\prime}=\varepsilon/(c\sqrt{N})$ which adds only a $\operatorname{poly}(\sqrt{N})$ factor), so we obtain the theorem. ∎

5. Proof of Proposition 1.3

In this section, we prove Proposition 1.3, showing that the density ratio for $MK^{-1/2}U_{NK}U_{NK}^{T}$ and $\mathbf{G}$ can be exponentially large when $K=\alpha M$ , $0<\alpha<1/2$ . We first consider $N=1$ , and prove

Lemma 5.1.

Let $K\geq 2$ and $M>K$ . The probability density function for $MK^{-1/2}U_{1K}U_{1K}^{T}$ over $\mathbb{C}$ is

(5.1)

\displaystyle f_{1}(z)=\frac{K(K-1)}{2\pi M^{2}B(K,M-K)}\int_{M^{-1}K^{1/2}|z|}^{1}(1-x)^{M-K-1}(x^{2}-M^{-2}K|z|^{2})^{(K-3)/2}\,dx,

where $B(x,y)=\frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)}$ denotes the beta function.

Proof.

The scalar quantity $U_{1K}U_{1K}^{T}$ is distributed as $u_{1}^{2}+\cdots+u_{K}^{2}$ for a random complex unit vector $u=(u_{1},\ldots,u_{M})\in\mathbb{S}_{\mathbb{C}}^{M-1}$ . For a vector $v$ , we will use notation $v_{a:b}$ to indicate the vector $(v_{a},v_{a+1},\ldots,v_{b})$ . Let $u=w/\|w\|_{2}$ for $w\sim\mathcal{CN}(0,I_{M})$ a standard complex Gaussian vector, so that

\displaystyle u_{1}^{2}+\cdots+u_{K}^{2}

\displaystyle=\frac{w_{1}^{2}+\cdots+w_{K}^{2}}{\|w_{1:K}\|_{2}^{2}}\frac{\|w_{1:K}\|_{2}^{2}}{\|w\|_{2}^{2}}=:SR,

where $S=\frac{w_{1}^{2}+\cdots+w_{K}^{2}}{\|w_{1:K}\|_{2}^{2}}\sim\sum_{j=1}^{K}s_{j}^{2}$ for $s=(s_{1},\ldots,s_{K})\sim\operatorname{Unif}(\mathbb{S}_{\mathbb{C}}^{K-1})$ , and $R:=\frac{\|w_{1:K}\|_{2}^{2}}{\|w\|_{2}^{2}}\sim\operatorname{Beta}(K,M-K)$ . Moreover, $S$ and $R$ are independent, since $(w_{1},\ldots,w_{K})/\|w_{1:K}\|_{2}$ is independent of $\|w_{1:K}\|_{2}$ , and also of $w_{K+1:M}$ . The density functions $f_{R}$ for $R\sim\operatorname{Beta}(K,M-K)$ on $[0,1]$ , and $f_{S}$ for $S$ on the unit disk, are

(5.2)

\displaystyle f_{R}(x)

\displaystyle=\frac{1}{B(K,M-K)}x^{K-1}(1-x)^{M-K-1},\;\text{ and }\;f_{S}(z)=\frac{K-1}{2\pi}(1-|z|^{2})^{(K-3)/2},

where $B(x,y)=\frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)}$ is the beta function. The density function $f_{S}$ was derived in [PM83, Eq. (3.12)]; alternatively, it can be derived by writing $S=\sum_{j=1}^{K}s_{j}^{2}$ , for $s=x+iy$ with $x,y\sim\mathcal{N}(0,I_{K})$ independent real Gaussian vectors, and then writing $|S|$ in terms of the eigenvalues of the real Wishart matrix $\begin{pmatrix}\|x\|_{2}^{2}&x\cdot y\\ x\cdot y&\|y\|_{2}^{2}\end{pmatrix}=VV^{T}$ for $V=\binom{-x-}{-y-}$ , since the density function for Wishart eigenvalues is well-known.

Since $S$ and $R$ are independent, the density function $f_{SR}$ for $SR$ on the unit disk is, e.g. by doing change of variables for $\mathbf{E}[h(SR)]$ for arbitrary $h$ ,

	$\displaystyle f_{SR}(z)$	$\displaystyle=\int_{\|z\|}^{1}dx\,f_{R}(x)f_{S}(z/x)\frac{1}{x^{2}}$
(5.3)			$\displaystyle=\frac{K-1}{2\pi B(K,M-K)}\int_{\|z\|}^{1}dx\,(1-x)^{M-K-1}(x^{2}-\|z\|^{2})^{(K-3)/2}.$

With scaling factors, the density function $f_{1}$ of $MK^{-1/2}U_{1K}U_{1K}^{T}$ is then given by (5.1). ∎

We can then give

Proof of Proposition 1.3.

We start with the case $N=1$ . The density function of $1\times 1$ $\mathbf{G}$ is $g_{1}(z)=\frac{1}{2\pi}e^{-|z|^{2}/2}$ , and the density of $MK^{-1/2}U_{1K}U_{1K}^{T}$ is given by (5.1) of Lemma 5.1.

For Laplace estimation, we consider $\frac{f_{1}(z)}{g_{1}(z)}=\frac{K(K-1)}{M^{2}B(K,M-K)}\int_{M^{-1}K^{1/2}|z|}^{1}e^{F_{z}(x)}\,dx$ , with

(5.4)

\displaystyle F_{z}(x)

\displaystyle=(M-K-1)\log(1-x)+\frac{K-3}{2}\log\left(x^{2}-M^{-2}K|z|^{2}\right)_{+}+\frac{|z|^{2}}{2}.

Solving $F_{z}^{\prime}(x)\equiv\frac{(K-3)x}{x^{2}-M^{-2}K|z|^{2}}-\frac{M-K-1}{1-x}=0$ and taking the positive root gives critical point

(5.5)

\displaystyle x_{z}

\displaystyle=\frac{K-3+\sqrt{(K-3)^{2}+4(M-4)(M-K-1)M^{-2}K|z|^{2}}}{2(M-4)}.

If $x_{z}$ is bounded away from the limits of integration as $M\to\infty$ , then one can check that Laplace’s method⁷⁷ 7 Since we only need a lower bound, we can actually skip Laplace’s method (to avoid having to check precise conditions for error bounds), and just integrate over an $O(M^{-1/2})$ neighborhood of $x_{*}$ , in which $F_{z_{*}}(x_{*})\geq F_{z_{*}}(x)-O(1)$ by Taylor’s theorem. gives

	$\displaystyle\frac{f_{1}(z)}{g_{1}(z)}$	$\displaystyle=\frac{K(K-1)}{M^{2}B(K,M-K)}\int_{M^{-1}K^{1/2}\|z\|}^{1}e^{F_{z}(x)}\,dx$
(5.6)			$\displaystyle=\frac{K(K-1)}{M^{2}B(K,M-K)}\sqrt{\frac{2\pi}{\|F_{z}^{\prime\prime}(x_{z})\|}}e^{F_{z}(x_{z})}\left(1+O(1/M)\right),$

with an implicit constant which may depend on $\alpha$ .

We can choose any $z$ with $|z|\leq MK^{-1/2}$ to try to make a large density ratio, so let’s solve for critical points of $F_{z}(x_{z})$ over $|z|$ . Writing $\phi(x,z):=F_{z}(x)$ , the multivariate chain rule gives $\frac{d}{d|z|}\phi(x_{z},z)=\partial_{x}\phi(x_{z},z)\frac{dx_{z}}{d|z|}+\partial_{z}\phi(x_{z},z)$ . Since $\partial_{x}\phi(x_{z},z)=F_{z}^{\prime}(x_{z})=0$ , we solve

\displaystyle(\partial_{z}F_{z})(x_{z})

\displaystyle\equiv\frac{K-3}{2}\frac{-2M^{-2}K|z|}{(x_{z}^{2}-M^{-2}K|z|^{2})}+|z|=0,

which has nonzero solutions

(5.7)

\displaystyle|z_{*}|^{2}:=M^{2}K^{-1}x_{z_{*}}^{2}-(K-3).

Let $x_{*}:=x_{z_{*}}$ , so $F^{\prime}_{z_{*}}(x_{*})=0$ . Plugging (5.7) into the equation $F^{\prime}_{z_{*}}(x_{*})=0$ implies

\displaystyle(K-3)x_{*}(1-x_{*})

\displaystyle=(M-K-1)M^{-2}K(K-3),

so $x_{*}(1-x_{*})=(M-K-1)M^{-2}K$ , and $x_{*}=\frac{1\pm\sqrt{1-4(M-K-1)M^{-2}K}}{2}$ . We will take the larger root and recall $\alpha<1/2$ , which gives

(5.8)		$\displaystyle x_{*}$	$\displaystyle=1-\alpha+\frac{\alpha}{(1-2\alpha)M}+O(1/M^{2}),$
(5.9)		$\displaystyle\|z_{*}\|^{2}$	$\displaystyle=\frac{M}{\alpha}(1-\alpha)^{2}-M\alpha+\frac{2(1-\alpha)}{1-2\alpha}+3+O(1/M).$

We also see that $M^{-1}K^{1/2}|z_{*}|=\sqrt{1-2\alpha}+O(1/M)$ , so for $\alpha>0$ there is a constant gap separating $x_{*}=1-\alpha+O(1/M)$ from the lower limit of integration in (5.6) as $M\to\infty$ . For this $z_{*}$ and $x_{*}$ ,

	$\displaystyle F_{z_{}}(x_{})$	$\displaystyle=(M-K-1)\log(1-x_{})+\frac{K-3}{2}\log(M^{-2}K(K-3))+\frac{1}{2}\|z_{}\|^{2}$
(5.10)			$\displaystyle=M\left[\log\alpha+\frac{1-2\alpha}{2\alpha}\right]-4\log\alpha+O(1/M).$

For the second derivative term in (5.6), evaluating the (negative of the) second derivative of $F_{z_{*}}$ at $x_{*}$ using (5.8) and (5.9) gives

	$\displaystyle-F_{z_{}}^{\prime\prime}(x_{})$	$\displaystyle=\frac{M-K-1}{(x_{}-1)^{2}}+\frac{(K-3)(x_{}^{2}+M^{-2}K\|z_{}\|^{2})}{(x_{}^{2}-M^{-2}K\|z_{*}\|^{2})^{2}}$
(5.11)			$\displaystyle=M\frac{2-3\alpha}{\alpha^{3}}+O(1).$

The remaining term to expand in (5.6) is the beta function $B(K,M-K)=\frac{\Gamma(K)\Gamma(M-K)}{\Gamma(M)}$ . Using Stirling’s formula $\log\Gamma(y)\sim(y-\frac{1}{2})\log y-y+\frac{1}{2}\log(2\pi)+O(1/y)$ , with $K,M-K\to\infty$ , we see

	$\displaystyle\log B$	$\displaystyle(K,M-K)$
		$\displaystyle=K\log\frac{K}{M}+(M-K)\log\frac{M-K}{M}+\frac{1}{2}\log(2\pi)+\frac{1}{2}\log\frac{M}{K(M-K)}+O(1/M)$
(5.12)			$\displaystyle=M\left[\alpha\log\alpha+(1-\alpha)\log(1-\alpha)\right]-\frac{1}{2}\log M-\frac{1}{2}\log(\alpha(1-\alpha))+\frac{1}{2}\log(2\pi)+O(1/M).$

In total, plugging (5.10), (5.11), and (5.12) into (5.6), we obtain

(5.13)

\displaystyle\log\frac{f_{1}(z_{*})}{g_{1}(z_{*})}

\displaystyle=M\left[(1-\alpha)\log\frac{\alpha}{1-\alpha}+\frac{1-2\alpha}{2\alpha}\right]+\frac{1}{2}\log(1-\alpha)-\frac{1}{2}\log(2-3\alpha)+O(1/M).

Thus

(5.14)

\displaystyle\sup_{z}\frac{f_{1}(z)}{g_{1}(z)}

\displaystyle\geq\sqrt{\frac{1-\alpha}{2-3\alpha}}\exp\left({M\left[(1-\alpha)\log\frac{\alpha}{1-\alpha}+\frac{1-2\alpha}{2\alpha}\right]}\right)(1-O(1/M)).

For $0<\alpha<1/2$ , we can check that

(5.15)

\displaystyle(1-\alpha)\log\frac{\alpha}{1-\alpha}+\frac{1-2\alpha}{2\alpha}>0,

for example by setting $y=\frac{1-\alpha}{\alpha}$ and noting the above is equivalent to $\frac{y^{2}-1}{2y}-\log y>0$ for $y>1$ , i.e. $\alpha<1/2$ , and that this holds by checking the derivative is $>0$ for $y>1$ . Then (5.14) gives

(5.16)

\displaystyle\sup_{z}\frac{f_{1}(z)}{g_{1}(z)}

\displaystyle\geq\Omega_{\alpha}(e^{c_{\alpha}M}),

with the implicit constant and $c_{\alpha}$ depending on $\alpha=K/M$ . Moreover, for $\alpha\in[\eta,1/2-\eta]$ , the constants can be taken uniform depending only on $\eta$ .

The $N=1$ case implies the higher dimensional cases since $f_{1}$ and $g_{1}$ are marginal densities of the higher dimensional densities $f$ and $g$ . If $\|f/g\|_{\infty}<\infty$ , then letting $d\hat{Z}$ denote integration over all variables except the top left entry $Z_{11}$ , we see

(5.17)

\displaystyle f_{1}(Z_{11})=\int d\hat{Z}\,\frac{f(Z)}{g(Z)}{g(Z)}

\displaystyle\leq\|f/g\|_{\infty}\int d\hat{Z}\,g(Z)=\|f/g\|_{\infty}g_{1}(Z_{11}).

Thus $\|f_{1}/g_{1}\|_{\infty}\leq\|f/g\|_{\infty}$ . ∎

Appendix A Proof of Theorem 1.1 with $N=o(\sqrt{K})$

In Section 2, we proved Proposition 2.1, which is a weaker version of Theorem 1.1. In this section, we prove the full Theorem 1.1. The general proof idea is still the same as for Proposition 2.1, in particular using the key Lemma 2.2. The difference is only in the later part of the proof, where we will do a more technical estimate to bound $d_{\mathrm{TV}}(Q\mathbf{G}Q^{T},\mathbf{G})$ .

Proposition 2.1 works for any $N=o(K^{1/3})$ , and the final bound (2.20) shows the proof also works for $N=o(K^{1/2})$ if $L=M-K=O(M^{1/2})$ . In this section, we improve the allowed size to $N=o(K^{1/2})$ for any $L$ . Since we already have the desired result $N=o(K^{1/2})$ when $L=O(M^{1/2})$ , and since $N=O(K^{1/2})=O(M^{1/2})$ , we will always consider $L\geq N$ in this section. It suffices to prove the required bound in Theorem 1.1 for $N^{2}\leq c_{0}K$ with $c_{0}$ chosen sufficiently small. For convenience, we will also take $K\geq 2N$ .

The proof starts in the same way as in Section 2, but the difference is we will bound the quantity $d_{\mathrm{TV}}(Q\mathbf{G}Q^{T},\mathbf{G})$ more sharply. Specifically, in this section we prove the stronger bound (compare to (2.20))

(A.1)

\displaystyle d_{\mathrm{TV}}(Q\mathbf{G}Q^{T},\mathbf{G})

\displaystyle\leq O\left(\frac{N}{\sqrt{K}}\right),

using an exact $\chi^{2}$ -divergence expression, and a standard Bakry–Émery concentration result [BE85, KM16]. Avoiding the inefficiency from Jensen’s inequality with the KL-divergence in (2.18) will allow us to obtain the error bound (A.1), at the cost of a more involved proof. Putting this new bound into (2.6) will then imply Theorem 1.1.

It will be useful to view symmetric $N\times N$ Gaussian matrices $\mathbf{G}$ and $q\mathbf{G}q^{T}$ directly as multivariate Gaussians in $d:=N(N+1)/2$ variables. To this end, let $\mathrm{Sym}_{N}:=\{X\in\mathbb{C}^{N\times N}:X^{T}=X\}\cong\mathbb{C}^{N(N+1)/2}$ , with the inner product $\langle X|Y\rangle=\frac{1}{2}\operatorname{Tr}(X^{\dagger}Y)$ . The usual orthonormal basis for $\mathrm{Sym}_{N}$ consists of $|i\rangle\langle j|+|j\rangle\langle i|$ for $1\leq i<j\leq N$ , and $\sqrt{2}|i\rangle\langle i|$ for $i=1,\ldots,N$ . For an invertible covariance operator $\Sigma:\mathrm{Sym}_{N}\to\mathrm{Sym}_{N}$ , the circular-symmetric complex Gaussian with covariance $\Sigma$ has density on $\mathrm{Sym}_{N}$ , with respect to the above orthonormal basis⁸⁸ 8 Note, due to the choice of inner product and resulting basis vector normalization, this differs by a constant factor from (2.7), which gives the density with respect to the upper triangular coordinates of the matrix., given by

(A.2)

\displaystyle f(Z)

\displaystyle=|\det\!{}_{\mathrm{Sym}_{N}}\Sigma|^{-1}\pi^{-d}e^{-\frac{1}{2}\operatorname{Tr}(Z^{\dagger}\Sigma^{-1}(Z))}.

We see, due to the orthonormal basis, that $\Sigma=I$ corresponds to the distribution of $\mathbf{G}\sim\mathcal{G}_{N}^{\mathrm{sym}}$ , whose entries have variance 2 on the diagonal and 1 off of the diagonal.

We can define a family of covariance operators $\Sigma_{R}$ on $\mathrm{Sym}_{N}$ via $\Sigma_{R}(X):=RXR^{T}$ , for $N\times N$ hermitian matrices $R$ . From a computation like (2.7), we see that for $q\geq 0$ , $q\mathbf{G}q^{T}$ has covariance operator $\Sigma_{q^{2}}$ . Also, when $R>0$ , then $\Sigma_{R}>0$ as well, since if $(r_{i},\varphi_{i})_{i}$ are eigenpairs of $R$ , then $(r_{i}r_{j},|\varphi_{i}\rangle\langle\bar{\varphi}_{j}|+|\varphi_{j}\rangle\langle\bar{\varphi}_{i}|)_{i\leq j}$ are eigenpairs of $\Sigma_{R}$ .

As in Section 2, let $P_{0}=I_{N}-E_{N}VV^{\dagger}E_{N}^{\dagger}$ , where $E_{N}=(I_{N}\;0_{M-N})$ and $V$ is the $M\times(M-K)$ matrix consisting of the last $M-K$ columns of the Haar random $U$ . Set $Q=(MK^{-1}P_{0})^{1/2}$ and $S=Q^{\dagger}Q=Q^{2}$ . For fixed $Q=q$ , $q\mathbf{G}q^{T}$ viewed as a Gaussian on $\mathrm{Sym}_{N}$ then has covariance $\Sigma_{s}$ , for $s=q^{\dagger}q=q^{2}$ . Write $s=q^{2}=:I+\delta$ as in Section 2, and let $\Sigma_{s}=I_{\mathrm{Sym}_{N}}+A$ with $A=A(\delta)=\Sigma_{s}-I_{\mathrm{Sym}_{N}}$ . As before, it will be useful to ensure $\|\delta\|$ is bounded away from 1. For random $Q=(MK^{-1}P_{0})^{1/2}$ and $\delta=Q^{2}-I$ , define the event

\displaystyle\mathcal{E}:=\{\|\delta\|\leq 1/4\},\quad\text{which has }\mathbf{P}[\mathcal{E}^{c}]=O\left(\frac{LN^{2}}{KM}\right),

by the same argument as (2.15). We will bound the distance between $\mu_{\mathcal{E}}=\mathcal{L}(Q\mathbf{G}Q^{T}|\mathcal{E})$ and $\gamma=\mathcal{L}(\mathbf{G})$ . Since $\mathbf{P}[\mathcal{E}^{c}]=O(\frac{LN^{2}}{KM})$ , this will also give us an adequate bound on the total variation distance between $\mu=\mathcal{L}(Q\mathbf{G}Q^{T})$ and $\mathcal{L}(\mathbf{G})$ .

It will be useful to use $\chi^{2}$ -divergence. The $\chi^{2}$ -divergence between probability measures $\mu$ and $\nu$ is $\chi^{2}(\mu||\nu):=\int\big(\frac{d\mu}{d\nu}-1\big)^{2}\,d\nu=\int\big(\frac{d\mu}{d\nu}\big)^{2}\,d\nu-1$ , and so $d_{\mathrm{TV}}(\mu,\nu)=\frac{1}{2}\int|\frac{d\mu}{d\nu}-1|\,d\nu\leq\frac{1}{2}\sqrt{\chi^{2}(\mu||\nu)}$ . The point is that for Gaussian mixture $\mu_{\mathcal{E}}$ and standard Gaussian $\gamma$ on $\mathrm{Sym}_{N}$ , we can evaluate an exact formula for the $\chi^{2}$ -divergence.

The density of $\mu_{\mathcal{E}}=\mathcal{L}(Q\mathbf{G}Q^{T}|\mathcal{E})$ is given as follows. Let $f_{A}$ be the density as in (A.2) for the complex Gaussian on $\mathrm{Sym}_{N}$ with covariance matrix $\Sigma=I_{\mathrm{Sym}_{N}}+A$ . Letting $A$ be distributed as the random variable $A(\delta)=\Sigma_{I+\delta}-I_{\mathrm{Sym}_{N}}$ conditioned on $\mathcal{E}$ , the density for $\mu_{\mathcal{E}}$ is then given by $\mathbf{E}_{A}f_{A}$ (for example, one can use Fubini–Tonelli). Then

	$\displaystyle 1+\chi^{2}(\mu_{\mathcal{E}}\|\|\gamma)=\int_{\mathrm{Sym}_{N}}\left(\frac{d\mu_{\mathcal{E}}}{d\gamma}\right)^{2}\,d\gamma$	$\displaystyle=\int_{\mathrm{Sym}_{N}}\frac{(\mathbf{E}_{A}f_{A}(Z))^{2}}{f_{0}^{2}(Z)}\,f_{0}(Z)\,dZ$
(A.3)			$\displaystyle=\int_{\mathrm{Sym}_{N}}\frac{\mathbf{E}_{A}f_{A}(Z)\mathbf{E}_{B}f_{B}(Z)}{f_{0}(Z)}\,dZ,$

where $B$ is an independent copy of $A$ . Since $f_{A},f_{B},f_{0}$ are Gaussian densities, we can evaluate, using (A.2),

	$\displaystyle 1+\chi^{2}(\mu_{\mathcal{E}}\|\|\gamma)$	$\displaystyle=\mathbf{E}_{A,B}\int_{\mathrm{Sym}_{N}}\pi^{-d}\det(I+A)^{-1}\det(I+B)^{-1}e^{-\frac{1}{2}\operatorname{Tr}(Z^{\dagger}[(I+A)^{-1}+(I+B)^{-1}-I](Z))}\,dZ$
(A.4)			$\displaystyle=\mathbf{E}_{A,B}\frac{\det(I+A)^{-1}\det(I+B)^{-1}}{\det((I+A)^{-1}+(I+B)^{-1}-I)}=\mathbf{E}_{A,B}\frac{1}{\det(I-AB)},$

the last equality by pulling out factors $\det(I+A)^{-1}$ and $\det(I+B)^{-1}$ on the left and right sides of the denominator, and provided all inverses are defined and $(I+A)^{-1}+(I+B)^{-1}-I>0$ , which we now verify.

To check $(I+A)^{-1}$ and $(I+B)^{-1}$ exist, we can estimate $\|A\|$ on $\mathcal{E}$ . Recall $S=I+\delta$ which is self-adjoint, and let $(\lambda_{j},\varphi_{j})_{j=1}^{N}$ be the eigenpairs of $\delta$ . Then one can check the (not necessarily normalized) eigenpairs of $\Sigma_{S}$ are $((1+\lambda_{i})(1+\lambda_{j}),|\varphi_{i}\rangle\langle\bar{\varphi}_{j}|+|\varphi_{j}\rangle\langle\bar{\varphi}_{i}|)_{i\leq j}$ . On $\mathcal{E}=\{\|\delta\|\leq 1/4\}$ , the absolute value of the eigenvalues of $A=\Sigma_{S}-I$ are thus

(A.5)

\displaystyle|(1+\lambda_{i})(1+\lambda_{j})-1|

\displaystyle=|\lambda_{i}+\lambda_{j}+\lambda_{i}\lambda_{j}|\leq\frac{9}{16}\quad\Longrightarrow\quad\|A\|\leq\frac{9}{16}<1.

Similarly, we can use the above bound to see $(I+A)^{-1}+(I+B)^{-1}-I\geq(\frac{16}{25}+\frac{16}{25}-1)I=\frac{7}{25}I>0$ . Thus we do have in total

(A.6)

\displaystyle 1+\chi^{2}(\mu_{\mathcal{E}}||\gamma)

\displaystyle=\mathbf{E}_{A,B}\det(I-AB)^{-1}.

The rest of the proof will be to bound $\mathbf{E}_{A,B}\det(I-AB)^{-1}$ . We will do this one expectation at a time. Letting $A=A(\delta)$ and taking a logarithm, define

(A.7)

\displaystyle F_{B}(\delta):=-\log\det(I-A(\delta)B),

so that $\mathbf{E}_{A}\det(I-AB)^{-1}=\mathbf{E}_{A}e^{F_{B}(\delta)}$ . Note that $\det(I-AB)$ is real and positive on $\mathcal{E}$ , using the factorization in (A.4) into ratios of determinants of positive operators on $\mathcal{E}$ .

We will use the Bakry–Émery criterion [BE85] to bound $\mathbf{E}_{A}e^{F_{B}(\delta)}$ via a concentration inequality. We use the version from [KM16, Theorem 2.1], which handles sets with not-necessarily-smooth boundaries.

Theorem A.1 (Bakry–Émery criterion [BE85, KM16]).

Suppose $\Omega$ is a convex subset of $\mathbb{R}^{m}$ (whose interior is thus geodesically convex, i.e. there is a distance-minimizing geodesic between any two points). Let $\mu$ be a probability measure of the form $d\mu(x)=e^{-V(x)}\,dx$ for $V$ smooth on the interior of $\Omega$ . Recall the entropy functional for nonnegative $h$ is $\operatorname{Ent}_{\mu}(h):=\int h\log h\,d\mu-\int h\,d\mu\cdot\log\int h\,d\mu$ . If $D^{2}V_{x}[H,H]\geq\lambda\|H\|_{\mathrm{hs}}^{2}$ for all $x\in\Omega$ , then for all $f\in C^{1}(\Omega)$ ,

(A.8)

\displaystyle\operatorname{Ent}_{\mu}(f^{2})

\displaystyle\leq\frac{2}{\lambda}\int|\nabla f|^{2}\,d\mu.

By the standard Herbst argument, the entropy bound with $f=e^{tF/2}$ implies for $L$ -Lipschitz $F$ (see e.g. [Wai19, §3.1.2]),

(A.9)

\displaystyle\log\mathbf{E}_{\mu}e^{t(F-\mathbf{E}F)}

\displaystyle\leq\frac{t^{2}L^{2}}{2\lambda},\quad t\in\mathbb{R}.

We apply this to $F=F_{B}$ and $V$ the negative log-density of $\delta$ restricted to $\Omega$ which will be a convex subset of $\mathcal{E}$ . We bound the Hessian of $V$ and the Lipschitz constant for $F_{B}$ as follows. For $V$ , recall $\delta=S-I=MK^{-1}P_{0}-I$ , where $P_{0}=(E_{N}A_{1})(E_{N}A_{1})^{\dagger}$ for $E_{N}A_{1}$ the top left $N\times K$ submatrix of the Haar random unitary $U$ . Then for $N\leq K,L$ , $P_{0}$ follows a complex matrix-variate beta distribution with density $f(P)\propto(\det P)^{K-N}\det(I-P)^{L-N}\mathbf{1}_{0<P<1}$ ; see the argument for the real orthogonal case in [Mec19, Lemma 2.12], which in the complex case gives the complex matrix-variate beta density in [DGGJ10, §2]. This gives the negative log-density of $\delta$ as

(A.10)

\displaystyle V(\delta)=-(K-N)\log\det(I+\delta)-(L-N)\log\det(L/K-\delta)+C^{\prime},

for some constant $C^{\prime}$ depending on $K,L,M,N$ , and with the domain $\mathbf{1}_{0<P<I}$ replaced by $\mathbf{1}_{-I<\delta<L/K}$ . Thus we take $\Omega:=\mathcal{E}\cap\{-I<\delta<L/K\}$ , which is convex.

The necessary matrix derivatives can be calculated using the formulas [PP12]

	$\displaystyle\partial_{t}\det A(t)$	$\displaystyle=\det A(t)\operatorname{Tr}(A(t)^{-1}\partial_{t}A(t)),$
	$\displaystyle\partial_{t}\operatorname{Tr}A(t)$	$\displaystyle=\operatorname{Tr}\partial_{t}A(t),$
	$\displaystyle\partial_{t}A^{-1}$	$\displaystyle=-A^{-1}(\partial_{t}A)A^{-1}.$

This gives

\displaystyle DV(\delta)[H]

\displaystyle=\partial_{t=0}V(\delta+tH)=-(K-N)\operatorname{Tr}((I+\delta)^{-1}H)+(L-N)\operatorname{Tr}((L/K-\delta)^{-1}H).

Taking another derivative using $\delta\mapsto\delta+tH$ , we get for $H$ hermitian,

	$\displaystyle D^{2}V(\delta)[H,H]$	$\displaystyle=(K-N)\operatorname{Tr}((I+\delta)^{-1}H(I+\delta)^{-1}H)+(L-N)\operatorname{Tr}((L/K-\delta)^{-1}H(L/K-\delta)^{-1}H)$
		$\displaystyle\geq(K-N)\operatorname{Tr}((I+\delta)^{-1}H(I+\delta)^{-1}H)$
(A.11)			$\displaystyle\geq(K-N)\left(\frac{4}{5}\right)^{2}\\|H\\|_{\mathrm{hs}}^{2}\geq cK\\|H\\|_{\mathrm{hs}}^{2},\quad\text{on $\mathcal{E}$,}$

and since we assume $K\geq 2N$ . For the Lipschitz constant for $F_{B}(\delta)=-\log\det(I_{\mathrm{Sym}_{N}}-A(\delta)B)$ , we similarly compute

(A.12)

\displaystyle DF_{B}(\delta)[H]

\displaystyle=\partial_{t=0}F_{B}(\delta+tH)=\operatorname{Tr}[(I-AB)^{-1}[DA(\delta)[H]]B].

Since we can also compute $DA(\delta)[H](X)=HX(1+\delta)^{T}+(1+\delta)XH^{T}$ , we can check that as an operator on $\mathrm{Sym}_{N}$ ,

(A.13)

\displaystyle\|DA(\delta)[H]\|_{\mathrm{hs}}\leq C\sqrt{N}\|H\|_{\mathrm{hs}},

for example by calculating the Hilbert–Schmidt norms in the usual basis for $\mathrm{Sym}_{N}$ . From (A.12), and $\|(I-AB)^{-1}\|\leq c$ and (A.5), we then get

(A.14)

\displaystyle\begin{aligned} \|\nabla F_{B}\|&=\sup_{H=H^{\dagger}:\|H\|_{\mathrm{hs}}=1}|DF_{B}(\delta)[H]|\leq C\sqrt{N}\|B\|_{\mathrm{hs}}.\end{aligned}

Combining with (A.11) and taking the convex set $\Omega=\mathcal{E}\cap\{-I<\delta<L/K\}$ , Theorem A.1 and (A.9) with $t=1$ thus give

(A.15)

\displaystyle\mathbf{E}_{A}e^{F_{B}}

\displaystyle\leq\exp\left[\mathbf{E}_{A}F_{B}+\frac{CN}{K}\|B\|_{\mathrm{hs}}^{2}\right].

To bound $\mathbf{E}_{A}F_{B}$ , we use the power series expansion $F_{B}(\delta)=-\log\det(I-A(\delta)B)=\sum_{j=1}^{\infty}\frac{1}{j}\operatorname{Tr}((AB)^{j})$ . For $j\geq 2$ ,

(A.16)

\displaystyle|\operatorname{Tr}((AB)^{j})|

\displaystyle\leq\|AB\|^{j-2}\|AB\|_{\mathrm{hs}}^{2}\leq\left(\frac{9}{16}\right)^{2(j-2)}\operatorname{Tr}(A^{2}B^{2}).

Thus $F_{B}(\delta)\leq\operatorname{Tr}(AB)+C\operatorname{Tr}(A^{2}B^{2})$ , and it will next be enough to bound $\mathbf{E}_{A}A$ and $\mathbf{E}_{A}A^{2}$ , which is done by the following lemma.

Lemma A.2 (expectation values).

Let $A=A(\delta)=\Sigma_{I+\delta}-I_{\mathrm{Sym}_{N}}$ , conditioned on the event $\mathcal{E}$ . Then there are scalars $a,b$ such that

(A.17)

\displaystyle\mathbf{E}A

\displaystyle=aI,\quad\mathbf{E}A^{2}=bI,\quad\text{with }a=O\bigg(\frac{\sqrt{N}}{K}\bigg),\quad b=O\left(\frac{N}{K}\right).

This lemma will be proved in Section A.1. Applying it to (A.15) and the discussion below it, we obtain

(A.18)

\displaystyle\mathbf{E}_{A}\det(I-AB)^{-1}

\displaystyle\leq e^{q(B)},\quad\text{ for }q(B)=a\operatorname{Tr}B+\left(b+\frac{CN}{K}\right)\operatorname{Tr}(B^{2}).

Now we apply Theorem A.1 with (A.9) again, this time to $\mathbf{E}_{B}e^{q(B)}$ . We calculate the gradient of $q(B)$ , letting $DB[H]$ denote $DB(\delta)[H]$ , which satisfies the bound (A.13),

(A.19)

\displaystyle Dq(B(\delta))[H]

\displaystyle=a\operatorname{Tr}(DB[H])+2(b+CN/K)\operatorname{Tr}[B(\delta)DB[H]].

Using (A.13), $b=O(\frac{N}{K})$ , $\|B\|_{\mathrm{hs}}\leq\sqrt{d}\|B\|=O(N)$ from (A.5), and $|\operatorname{Tr}_{\mathrm{Sym}_{N}}Y|\leq cN\|Y\|_{\mathrm{hs}}$ and $|\operatorname{Tr}XY|\leq\|X\|_{\mathrm{hs}}\|Y\|_{\mathrm{hs}}$ , we see on $\mathcal{E}$ ,

(A.20)

\displaystyle\|\nabla q(B)\|

\displaystyle\leq O\left(|a|N^{3/2}+\frac{N^{5/2}}{K}\right).

Then Theorem A.1 with (A.9), and $a=O(\frac{\sqrt{N}}{K})$ , give

(A.21)

\displaystyle\mathbf{E}_{B}e^{q(B)}

\displaystyle\leq\exp\left[\mathbf{E}_{B}q(B)+\frac{C}{K}\left(\frac{N^{4}}{K^{2}}+\frac{N^{5}}{K^{2}}\right)\right].

Using Lemma A.2 for $B$ , along with the bounds on $a=O(\frac{\sqrt{N}}{K})$ and $b=O(\frac{N}{K})$ , then gives

(A.22)

\displaystyle\mathbf{E}_{B}q(B)

\displaystyle\leq da^{2}+d\frac{C^{\prime}N}{K}b\leq O\left(\frac{N^{3}}{K^{2}}+\frac{N^{4}}{K^{2}}\right)=O\left(\frac{N^{4}}{K^{2}}\right),

since $N^{2}\leq c_{0}K$ by assumption. Thus

(A.23)

\displaystyle\mathbf{E}_{A,B}\det(I-AB)^{-1}

\displaystyle\leq\mathbf{E}_{B}e^{q(B)}\leq\exp\left(O\bigg(\frac{N^{4}}{K^{2}}\bigg)\right).

Finally, returning to (A.6), we get for $N^{2}\leq c_{0}K$ ,

(A.24)

\displaystyle\chi^{2}(\mu_{\mathcal{E}}||\gamma)

\displaystyle\leq e^{O\left(\frac{N^{4}}{K^{2}}\right)}-1\leq O\left(\frac{N^{4}}{K^{2}}\right).

Since $d_{\mathrm{TV}}(\mu_{\mathcal{E}},\gamma)\leq\frac{1}{2}\sqrt{\chi^{2}(\mu_{\mathcal{E}}||\gamma)}$ , and $d_{\mathrm{TV}}(\mu_{\mathcal{E}},\mu)\leq\mathbf{P}[\mathcal{E}^{c}]\leq O(\frac{LN^{2}}{KM})$ for $\mu=\mathcal{L}(Q\mathbf{G}Q^{T})$ (for example, use a coupling $Z=Q\mathbf{G}Q^{T}\sim\mu$ , $Y=Z$ on $\mathcal{E}(Q)$ , $Y\sim\mu_{\mathcal{E}}$ on $\mathcal{E}(Q)^{c}$ ), we obtain

(A.25)

\displaystyle d_{\mathrm{TV}}(Q\mathbf{G}Q^{T},\mathbf{G})

\displaystyle\leq O\left(\frac{N}{\sqrt{K}}\right).

Inserting this in (2.6) gives Theorem 1.1. ∎

A.1. Proof of Lemma A.2

In this section, we prove Lemma A.2. It will be an application of the following result, followed by some trace bounds similar to those in Section 2.

Lemma A.3.

Let $R$ be a random hermitian matrix with law invariant under $R\mapsto VRV^{\dagger}$ for any $V\in\mathrm{U}(N)$ . Then $\mathbf{E}\Sigma_{R}=a_{R}I$ for some $a_{R}\in\mathbb{R}$ .

Proof of Lemma A.3.

It will be enough to show that $\mathbf{E}\Sigma_{R}(|x\rangle\langle x|)=a|x\rangle\langle x|$ for every $x\in\mathbb{R}^{N}$ , since the real basis elements $|e_{i}\rangle\langle e_{j}|+|e_{j}\rangle\langle e_{i}|$ for $\mathrm{Sym}_{N}$ can be expressed as $|e_{i}+e_{j}\rangle\langle e_{i}+e_{j}|-|e_{i}\rangle\langle e_{i}|-|e_{j}\rangle\langle e_{j}|$ . We start with $|x\rangle\langle x|=|e_{1}\rangle\langle e_{1}|=\begin{pmatrix}1&\mathbf{0}\\ \mathbf{0}&\mathbf{0}_{N-1}\end{pmatrix}$ . We see that $\mathbf{E}\Sigma_{R}(|e_{1}\rangle\langle e_{1}|)$ is invariant under the unitary congruence $(\cdot)\mapsto D(\cdot)D^{T}$ by the block matrix $D:=\operatorname{diag}(1,W)$ for any $W\in\mathrm{U}(N-1)$ , since noting that $D|e_{1}\rangle\langle e_{1}|D^{T}=|e_{1}\rangle\langle e_{1}|$ , then

	$\displaystyle DR\|e_{1}\rangle\langle e_{1}\|R^{T}D^{T}$	$\displaystyle=DRD^{\dagger}D\|e_{1}\rangle\langle e_{1}\|D^{T}(D^{T})^{\dagger}R^{T}D^{T}$
(A.26)			$\displaystyle\overset{d}{=}R\|e_{1}\rangle\langle e_{1}\|R^{T}.$

Write $\mathbf{E}\Sigma_{R}(|e_{1}\rangle\langle e_{1}|)\equiv\mathbf{E}R|e_{1}\rangle\langle e_{1}|R^{T}=\begin{pmatrix}a&b\\ b^{T}&C\end{pmatrix}$ for $a\in\mathbb{R}$ and $C$ an $(N-1)\times(N-1)$ complex symmetric matrix. Then taking $W=e^{it}I_{N-1}$ for $t\in\mathbb{R}$ in (A.26) shows $b=0$ and $C=0$ , so $\mathbf{E}\Sigma_{R}(|e_{1}\rangle\langle e_{1}|)=a|e_{1}\rangle\langle e_{1}|$ . For general unit $x\in\mathbb{R}^{N}$ , we can then simply let $V\in\mathrm{O}(N)$ be an orthogonal matrix such that $V|e_{1}\rangle=|x\rangle$ ; then $V^{\dagger}=V^{T}$ and

(A.27)

\displaystyle\Sigma_{R}(|x\rangle\langle x|)=R|x\rangle\langle x|R^{T}

\displaystyle=RV|e_{1}\rangle\langle e_{1}|V^{T}R^{T}\overset{d}{=}VR|e_{1}\rangle\langle e_{1}|R^{T}V^{T},

and so $\mathbf{E}\Sigma_{R}(|x\rangle\langle x|)=Va|e_{1}\rangle\langle e_{1}|V^{T}=a|x\rangle\langle x|$ , as desired. ∎

We apply this to estimate $\mathbf{E}_{\mathcal{E}}A$ and $\mathbf{E}_{\mathcal{E}}A^{2}$ in Lemma A.2. Since $S=1+\delta=MK^{-1}P_{0}=MK^{-1}(I_{N}-E_{N}VV^{\dagger}E_{N}^{\dagger})$ for $V$ the last $M-K$ columns of an $M\times M$ Haar random matrix, one can check that $S\overset{d}{=}\xi_{N}S\xi_{N}^{\dagger}$ for any $\xi_{N}\in\mathrm{U}(N)$ . (For example, one can use $\xi_{N}E_{N}=[\xi_{N}\;0_{M-N}]=E_{N}(\xi_{N}\oplus I_{M-N})$ , following by invariance of $V$ under multiplication by $(\xi_{N}\oplus I_{M-N})\in\mathrm{U}(M)$ .) Moreover, since conjugation preserves eigenvalues and $S=I+\delta$ , the distributions conditioned on $\mathcal{E}$ are also the same. From the lemma, then $\mathbf{E}_{\mathcal{E}}A=\mathbf{E}_{\mathcal{E}}\Sigma_{S}-I=aI$ for some $a$ , and $\mathbf{E}_{\mathcal{E}}A^{2}=\mathbf{E}_{\mathcal{E}}\Sigma_{S^{2}}-2\mathbf{E}_{\mathcal{E}}\Sigma_{S}+I=bI$ for some $b$ .

•

For $\mathbf{E}_{\mathcal{E}}A^{2}=bI$ , in terms of the eigenvalues $\lambda_{j}$ for $\delta$ (see (A.5)), we have $\operatorname{Tr}_{\mathrm{Sym}_{N}}(A^{2})=\sum_{i\leq j}(\lambda_{i}+\lambda_{j}+\lambda_{i}\lambda_{j})^{2}\leq\sum_{i\leq j}C(\lambda_{i}^{2}+\lambda_{j}^{2})$ . Thus

	$\displaystyle b=\frac{1}{\dim\mathrm{Sym}_{N}}\mathbf{E}_{\mathcal{E}}\operatorname{Tr}_{\mathrm{Sym}_{N}}(A^{2})$	$\displaystyle\leq\frac{2}{N(N+1)}\mathbf{E}_{\mathcal{E}}\sum_{i\leq j}C(\lambda_{i}^{2}+\lambda_{j}^{2})$
(A.28)			$\displaystyle\leq\frac{C^{\prime}}{N+1}\mathbf{E}_{\mathcal{E}}\operatorname{Tr}\delta^{2}=O\left(\frac{LN}{KM}\right),$

where we used (2.14) and that $\mathbf{P}[\mathcal{E}]\geq 1-c$ to conclude for any $X\geq 0$ , $\mathbf{E}[X|\mathcal{E}]=\mathbf{E}[X\mathbf{1}_{\mathcal{E}}]/\mathbf{P}[\mathcal{E}]\leq C\mathbf{E}[X]$ .

•

For $\mathbf{E}_{\mathcal{E}}A=aI$ , we have $a=\frac{2}{N(N+1)}\mathbf{E}_{\mathcal{E}}\operatorname{Tr}_{\mathrm{Sym}_{N}}A$ . Letting $(s_{i})_{i}$ be the eigenvalues of $S$ , then the eigenvalues of $\Sigma_{S}$ are $(s_{i}s_{j})_{i\leq j}$ (see the discussion before (A.5)), and so

(A.29)

\displaystyle\operatorname{Tr}_{\mathrm{Sym}_{N}}\Sigma_{S}

\displaystyle=\sum_{i\leq j}s_{i}s_{j}=\frac{1}{2}[(\operatorname{Tr}S)^{2}+\operatorname{Tr}(S^{2})].

Since $A=\Sigma_{S}-I$ and $S=I+\delta$ , we then have

(A.30)

\displaystyle\operatorname{Tr}_{\mathrm{Sym}_{N}}A=\frac{1}{2}[(\operatorname{Tr}\delta)^{2}+\operatorname{Tr}(\delta^{2})]+(N+1)\operatorname{Tr}\delta.

Using $\mathbf{E}\operatorname{Tr}\delta=0$ , we can bound $\mathbf{E}_{\mathcal{E}}\operatorname{Tr}\delta$ as

\displaystyle|\mathbf{E}_{\mathcal{E}}\operatorname{Tr}\delta|

\displaystyle=\frac{|\mathbf{E}\operatorname{Tr}\delta\mathbf{1}_{\mathcal{E}^{c}}|}{\mathbf{P}[\mathcal{E}]}\leq C\sqrt{\mathbf{E}(\operatorname{Tr}\delta)^{2}}\sqrt{\mathbf{P}[\mathcal{E}^{c}]}.

Dividing (A.30) by $\dim\mathrm{Sym}_{N}=\frac{N(N+1)}{2}$ and using (2.13), (2.14) to bound the unconditioned expectation values of $(\operatorname{Tr}\delta)^{2}$ and $\operatorname{Tr}(\delta^{2})$ terms, and (2.15) for $\mathbf{P}[\mathcal{E}^{c}]\leq O(\frac{LN^{2}}{KM})$ , we get

(A.31)

\displaystyle|a|\leq O\bigg(\frac{L\sqrt{N}}{KM}\bigg)=O\bigg(\frac{\sqrt{N}}{K}\bigg).

This proves Lemma A.2. ∎

Acknowledgments. This project used GPT-5.5 Thinking and Pro and GPT-5.6 Sol for coming up with proof ideas and methods, as well as for general checking and proofreading. The paper was written by the authors and all results and proofs were checked and validated by the authors, who are fully responsible for the final content. We thank Joseph Iosue and Yu-Xin Wang for useful discussions. L.S. and A.V.G. acknowledge support from the U.S. Department of Energy, Office of Science, Accelerated Research in Quantum Computing, Fundamental Algorithmic Research toward Quantum Utility (FAR-Qu). L.S. and A.V.G. were also supported in part by ARL (W911NF-24-2-0107), ONR MURI, and NSF QLCI (award No. OMA-2120757). V.G. was supported by the US Army Research Office under Grant Number W911NF-23-1-0241.

References

[AA13] S. Aaronson and A. Arkhipov, The computational complexity of linear optics, Theory of Computing 9 (2013), no. 4, 143–252.
[Aar03] S. Aaronson, Algorithms for boolean function query properties, SIAM Journal on Computing 32 (2003), no. 5, 1140–1157.
[BBD⁺26] A. Bouland, D. Brod, I. Datta, B. Fefferman, D. Grier, F. Hernández, and M. Oszmaniec, Complexity-theoretic foundations of BosonSampling with a linear number of modes, Phys. Rev. X 16 (2026), 021059.
[BDER16] S. Bubeck, J. Ding, R. Eldan, and M. Z. Rácz, Testing for high-dimensional geometry in random graphs, Random Structures Algorithms 49 (2016), no. 3, 503–532.
[BE85] D. Bakry and M. Émery, Diffusions hypercontractives, Séminaire de probabilités, XIX, 1983/84, Lecture Notes in Math., vol. 1123, Springer, Berlin, 1985, pp. 177–206.
[BGP00] M. Bellare, O. Goldreich, and E. Petrank, Uniform generation of NP-witnesses using an NP-oracle, Information and Computation 163 (2000), no. 2, 510–526.
[Col03] B. Collins, Moments and cumulants of polynomial random variables on unitary groups, the Itzykson-Zuber integral, and free probability, Int. Math. Res. Not. (2003), no. 17, 953–982.
[com] Complexity zoo, https://complexityzoo.net/Complexity˙Zoo, Accessed: 2026-06.
[CS06] B. Collins and P. Śniady, Integration with respect to the Haar measure on unitary, orthogonal and symplectic group, Comm. Math. Phys. 264 (2006), no. 3, 773–795.
[DGGJ10] J. A. Díaz-García and R. Gutiérrez Jáimez, Complex bimatrix variate generalised beta distributions, Linear Algebra Appl. 432 (2010), no. 2-3, 571–582.
[DGL⁺23] Y.-H. Deng, Y.-C. Gu, H.-L. Liu, S.-Q. Gong, H. Su, Z.-J. Zhang, H.-Y. Tang, M.-H. Jia, J.-M. Xu, M.-C. Chen, J. Qin, L.-C. Peng, J. Yan, Y. Hu, J. Huang, H. Li, Y. Li, Y. Chen, X. Jiang, L. Gan, G. Yang, L. You, L. Li, H.-S. Zhong, H. Wang, N.-L. Liu, J. J. Renema, C.-Y. Lu, and J.-W. Pan, Gaussian boson sampling with pseudo-photon-number-resolving detectors and quantum computational advantage, Phys. Rev. Lett. 131 (2023), 150601.
[DMV⁺22] A. Deshpande, A. Mehta, T. Vincent, N. Quesada, M. Hinsche, M. Ioannou, L. Madsen, J. Lavoie, H. Qi, J. Eisert, D. Hangleiter, B. Fefferman, and I. Dhand, Quantum computational advantage via high-dimensional Gaussian boson sampling, Science Advances 8 (2022), no. 1, eabi7894.
[EID⁺25a] A. Ehrenberg, J. T. Iosue, A. Deshpande, D. Hangleiter, and A. V. Gorshkov, Second moment of hafnians in Gaussian boson sampling, Phys. Rev. A 111 (2025), 042412.
[EID⁺25b] A. Ehrenberg, J. T. Iosue, A. Deshpande, D. Hangleiter, and A. V. Gorshkov, Transition of anticoncentration in Gaussian boson sampling, Phys. Rev. Lett. 134 (2025), 140601.
[GBA⁺22] D. Grier, D. J. Brod, J. M. Arrazola, M. B. de Andrade Alonso, and N. Quesada, The complexity of bipartite Gaussian boson sampling, Quantum 6 (2022), 863.
[HHT97] Y. Han, L. A. Hemaspaandra, and T. Thierauf, Threshold computation and cryptographic security, SIAM J. Comput. 26 (1997), no. 1, 59–78.
[HKS⁺17] C. S. Hamilton, R. Kruse, L. Sansoni, S. Barkhofen, C. Silberhorn, and I. Jex, Gaussian boson sampling, Phys. Rev. Lett. 119 (2017), 170501.
[JL15] T. Jiang and D. Li, Approximation of rectangular beta-Laguerre ensembles and large deviations, J. Theoret. Probab. 28 (2015), no. 3, 804–847.
[KHS⁺19] R. Kruse, C. S. Hamilton, L. Sansoni, S. Barkhofen, C. Silberhorn, and I. Jex, Detailed study of Gaussian boson sampling, Phys. Rev. A 100 (2019), 032326.
[KM16] A. V. Kolesnikov and E. Milman, Riemannian metrics on convex sets with applications to Poincaré and log-Sobolev inequalities, Calc. Var. Partial Differential Equations 55 (2016), no. 4, Art. 77, 36.
[LSD⁺26] H.-L. Liu, H. Su, Y.-H. Deng, S.-Q. Gong, Y.-C. Gu, H.-Y. Tang, M.-H. Jia, Q. Wei, Y.-K. Song, D.-Z. Wang, M.-Y. Zheng, F.-X. Chen, L.-B. Li, S.-Y. Ren, X.-Z. Zhu, M.-H. Wang, Y.-J. Chen, Y.-F. Liu, L.-S. Song, P.-Y. Yang, J.-S. Chen, H. An, L. Zhang, L. Gan, G.-w. Yang, J.-M. Xu, Y.-M. He, H. Wang, H.-S. Zhong, M.-C. Chen, X. Jiang, L. Li, N.-L. Liu, X.-L. Su, Q. Zhang, C.-Y. Lu, and J.-W. Pan, Gaussian boson sampling with 1,024 squeezed states in 8,176 modes, Nature 653 (2026), no. 8115, 687–692.
[Mec19] E. S. Meckes, The random matrix theory of the classical compact groups, Cambridge Tracts in Mathematics, vol. 218, Cambridge University Press, Cambridge, 2019.
[MLA⁺22] L. S. Madsen, F. Laudenbach, M. F. Askarani, F. Rortais, T. Vincent, J. F. Bulmer, F. M. Miatto, L. Neuhaus, L. G. Helt, M. J. Collins, et al., Quantum computational advantage with a programmable photonic processor, Nature 606 (2022), no. 7912, 75–81.
[PM83] P. Pereyra and P. A. Mello, Marginal distribution of the $S$ -matrix elements for Dyson’s measure and some applications, J. Phys. A 16 (1983), no. 2, 237–254.
[PP12] K. B. Petersen and M. S. Pedersen, The matrix cookbook, 2012.
[RR19] M. Z. Rácz and J. Richey, A smooth transition from Wishart to GOE, J. Theoret. Probab. 32 (2019), no. 2, 898–906.
[SMG26] L. Shou, S. H. Miller, and V. Galitski, Proof of hiding conjecture in Gaussian boson sampling, to appear in Quantum (2026), arXiv:2508.00983.
[Sto85] L. Stockmeyer, On approximation algorithms for #P, SIAM Journal on Computing 14 (1985), no. 4, 849–861.
[Tod91] S. Toda, PP is as hard as the polynomial-time hierarchy, SIAM J. Comput. 20 (1991), no. 5, 865–877.
[Val79] L. Valiant, The complexity of computing the permanent, Theoretical Computer Science 8 (1979), no. 2, 189–201.
[vH16] R. van Handel, Probability in High Dimension, APC 550 Lecture Notes, 2016.
[Wai19] M. J. Wainwright, High-dimensional statistics, Cambridge Series in Statistical and Probabilistic Mathematics, vol. 48, Cambridge University Press, Cambridge, 2019.
[ZDQ⁺21] H.-S. Zhong, Y.-H. Deng, J. Qin, H. Wang, M.-C. Chen, L.-C. Peng, Y.-H. Luo, D. Wu, S.-Q. Gong, H. Su, Y. Hu, P. Hu, X.-Y. Yang, W.-J. Zhang, H. Li, Y. Li, X. Jiang, L. Gan, G. Yang, L. You, Z. Wang, L. Li, N.-L. Liu, J. J. Renema, C.-Y. Lu, and J.-W. Pan, Phase-programmable Gaussian boson sampling using stimulated squeezed light, Phys. Rev. Lett. 127 (2021), 180502.
[ZWD⁺20] H.-S. Zhong, H. Wang, Y.-H. Deng, M.-C. Chen, L.-C. Peng, Y.-H. Luo, J. Qin, D. Wu, X. Ding, Y. Hu, P. Hu, X.-Y. Yang, W.-J. Zhang, H. Li, Y. Li, X. Jiang, L. Gan, G. Yang, L. You, Z. Wang, L. Li, N.-L. Liu, C.-Y. Lu, and J.-W. Pan, Quantum computational advantage using photons, Science 370 (2020), no. 6523, 1460–1463.

	$\displaystyle\left\|\|\operatorname{Haf}X\|^{2}-\|\operatorname{Haf}Y\|^{2}\right\|$	$\displaystyle=\Bigg\|\sum_{\pi,\pi^{\prime}\in\mathcal{P}_{2}(N)}\left(\prod_{\{i,j\}\in\pi}X_{ij}\prod_{\{i^{\prime},j^{\prime}\}\in\pi^{\prime}}\bar{X}_{i^{\prime}j^{\prime}}-\prod_{\{i,j\}\in\pi}Y_{ij}\prod_{\{i^{\prime},j^{\prime}\}\in\pi^{\prime}}\bar{Y}_{i^{\prime}j^{\prime}}\right)\Bigg\|$
		$\displaystyle\leq\sum_{\pi,\pi^{\prime}\in\mathcal{P}_{2}(N)}NB^{N-1}\\|X-Y\\|_{\infty}$
(4.14)			$\displaystyle\leq(N-1)!!^{2}NB^{N-1}\\|X-Y\\|_{\infty}.$

	$\displaystyle\mathbf{P}\left[\|\tilde{q}_{S^{}}-p_{S^{}}\|>\varepsilon\frac{N!}{M^{N}}\right]$	$\displaystyle\leq\mathbf{P}\left[\|\tilde{q}_{S^{}}-q_{S^{}}\|>\frac{\varepsilon}{2}\frac{N!}{M^{N}}\right]+\mathbf{P}\left[\|q_{S^{}}-p_{S^{}}\|>\frac{\varepsilon}{2}\frac{N!}{M^{N}}\right]$
		$\displaystyle\leq\mathbf{P}\left[\|\tilde{q}_{S^{}}-q_{S^{}}\|>\alpha q_{S^{}}\right]+\mathbf{P}\left[q_{S^{}}>2j\frac{N!}{M^{N}}\right]+O(\delta)$
(4.24)			$\displaystyle\leq\frac{1}{2^{M}}+\frac{\delta}{4}+O(\delta).$

	$\displaystyle DR\|e_{1}\rangle\langle e_{1}\|R^{T}D^{T}$	$\displaystyle=DRD^{\dagger}D\|e_{1}\rangle\langle e_{1}\|D^{T}(D^{T})^{\dagger}R^{T}D^{T}$
(A.26)			$\displaystyle\overset{d}{=}R\|e_{1}\rangle\langle e_{1}\|R^{T}.$

Proof of the hiding conjecture for Gaussian boson sampling with an arbitrary number of squeezed input modes

Abstract.

1. Introduction

1.1. Preliminaries

Conjecture 1 (hiding in Gaussian boson sampling).

1.2. Main results

Theorem 1.1 (hiding in Gaussian boson sampling).

Remark 1.1.

Theorem 1.2 (G​GTGG^{T} vs 𝒢Nsym\mathcal{G}_{N}^{\mathrm{sym}}).

Remark 1.2.

Proposition 1.3 (no density ratio bound).

Remark 1.3.

Problem 1 (|GHE|±2|\mathrm{GHE}|_{\pm}^{2}).

Theorem 1.4 (main hardness result).

Conjecture 2.

1.3. Outline

2. Proof of N=o⁡(K1/3)N=o(K^{1/3}) TVD hiding

Proposition 2.1 (hiding for N=o⁡(K1/3)N=o(K^{1/3})).

Lemma 2.2.

Proof.

Proof of Proposition 2.1.

3. Proof of Theorem 1.2

4. Hardness argument

Assumption 1.

Lemma 4.1 (approximate instance generation).

Proof of Lemma 4.1.

Lemma 4.2 (boundary estimate).

Proof.

Proof of Theorem 1.4.

5. Proof of Proposition 1.3

Lemma 5.1.

Proof.

Proof of Proposition 1.3.

Appendix A Proof of Theorem 1.1 with N=o⁡(K)N=o(\sqrt{K})

Theorem A.1 (Bakry–Émery criterion [BE85, KM16]).

Lemma A.2 (expectation values).

A.1. Proof of Lemma A.2

Lemma A.3.

Proof of Lemma A.3.

References

Theorem 1.2 ( $GG^{T}$ vs $\mathcal{G}_{N}^{\mathrm{sym}}$ ).

Problem 1 ( $|\mathrm{GHE}|_{\pm}^{2}$ ).

2. Proof of $N=o(K^{1/3})$ TVD hiding

Proposition 2.1 (hiding for $N=o(K^{1/3})$ ).

Appendix A Proof of Theorem 1.1 with $N=o(\sqrt{K})$