# Classical and Quantum Algorithms for Tensor Principal Component Analysis

Matthew B. Hastings  Station Q, Microsoft Research, Santa Barbara, CA
93106-6105, USA  Microsoft Quantum and Microsoft Research, Redmond, WA 98052,
USA

###### Abstract

We present classical and quantum algorithms based on spectral methods for a
problem in tensor principal component analysis. The quantum algorithm achieves
a quartic speedup while using exponentially smaller space than the fastest
classical spectral algorithm, and a super-polynomial speedup over classical
algorithms that use only polynomial space. The classical algorithms that we
present are related to, but slightly different from those presented recently
in Ref. [1]. In particular, we have an improved threshold for recovery and the
algorithms we present work for both even and odd order tensors. These results
suggest that large-scale inference problems are a promising future application
for quantum computers.

##  1 Introduction

Principal component analysis is a fundamental technique that finds
applications in reducing the dimensionality of data and denoising. While an
optimal choice of principal components for a matrix can be computed
efficiently using linear algebra, the corresponding problem for tensors is
much less well-understood. Ref. [2] introduced a simple statistical model for
tensor principal component analysis, termed the “spiked tensor” problem, and
this paper has lead to a large amount of follow-up research. The model
consists of (see below for more precise definitions) randomly choosing some
unknown “signal vector” $v_{\rm sig}\in{\mathbb{R}}^{N}$; then, the $p$-th
order tensor

$T_{0}=\lambda v_{\rm sig}^{\otimes p}+G$ (1)

is formed, where $G\in({\mathbb{R}^{N}})^{\otimes p}$ is noise chosen from
some random distribution and where $\lambda$ is some scalar representing a
signal-to-noise ratio. One task called recovery is to infer $v_{\rm sig}$ (to
some accuracy) given $T_{0}$. A simpler task called detection is to
distinguish the case $\lambda=0$ from $\lambda=\overline{\lambda}$ for some
$\overline{\lambda}>0$, again just given $T_{0}$.

Ref. [2] presented a variety of algorithms for this problem. Following the
normalization of Ref. [1], the entries of $G$ are chosen independently from a
Gaussian distribution of zero mean and unit variance, with $|v_{\rm
sig}|=\sqrt{N}$. Then, information theoretically, it is possible to recover
for $\lambda$ much larger than $N^{(1-p)/2}$ [2, 3]. However, no polynomial
time algorithm is known that achieves this performance. Rather, the two best
known algorithms are spectral and sum-of-squares. Spectral algorithms were
first suggested in Ref. [2]. There, a matrix is formed from $T_{0}$ (if $p$ is
even, the matrix is $N^{p/2}$-by-$N^{p/2}$, with its entries given by entries
of $T_{0}$) and the leading eigenvector of the matrix is used to determine
$v_{\rm sig}$. For even $p$, this method works for $\lambda$ much larger than
$N^{-p/4}$, and a variant of it is conjectured to perform similarly for odd
$p$. Methods based on the sum-of-squares also perform similarly to the
spectral method. The sum-of-squares method [4, 5] for this problem gives rise
to a sequence of algorithms [6, 7], in which one can recover at $\lambda$
smaller than $N^{-p/4}$ at the cost of runtime and space increasing
exponentially in ${\rm polylog}(N)N^{-p/4}/\lambda$. In Ref. [1], a sequence
of spectral algorithms with similar performance was shown.

In this paper, we present another spectral algorithm for the problem. Our
spectral algorithm for even $p$ is closely related to that of Ref. [1] which
we became aware of while this paper was in preparation, and we use the
normalization in that paper. However, we make several changes. Our main
technical results are the following. First, we prove an improved threshold for
recovery for even $p$ using our algorithm; the improvement is by a constant
factor and relies on a randomized method of recovering. Second, we provide a
different algorithm for odd $p$ with provable guarantees on while no
guarantees were given in Ref. [1] for odd $p$. For both even and odd $p$, we
have provable bounds on recovery for $\lambda$ of order $N^{-p/4}$ (without
any polylogarithmic factors) and we have a sequence of algorithms similar to
that above for $\lambda$ small compared to $N^{-p/4}$. Third, we give a
quantum algorithm for our spectral method, achieving a quartic speedup and
exponential reduction in space. This quantum algorithm involves two main
ideas. The first uses phase estimation and amplitude amplification to obtain a
quadratic speedup in computing the largest eigenvector. The second idea uses a
chosen input state to obtain a further quadratic speedup, giving the overall
quartic speedup.

We emphasize that the quantum speedup is quartic compared to classical
spectral algorithms presented here and in previous work. We are not able to
make an accurate comparison of the runtime to sum-of-squares methods. In part,
given that the runtime of all of these algorithms increases exponentially in
$\lambda^{-1}$, a change in prefactors in some estimates for threshold can
give rise to polynomial changes in runtime. We expect that many of these
estimates of thresholds are not tight (indeed, we expect that they are off by
a polylogarithmic factor), and so either improved analytic methods or
numerical simulations are needed to give an accurate comparison.

At a more heuristic level, we present a rather different motivation for our
spectral algorithm compared to Ref. [1]. Rather than being motivated by the
so-called Kikuchi free energy, we instead are motivated by mean-field
approximations to quantum many-body systems. We consider a system of a some
number ${n_{bos}}$ of qudits, each of dimension $N$, and use the tensor
$T_{0}$ to construct a quantum Hamiltonian on these qudits. Increasing
${n_{bos}}$ gives rise to a similar sequence of algorithms as above, with
increased runtime but improved performance: the required ${n_{bos}}$ increases
polynomially in $\lambda^{-1}$ as ${\rm
polylog}(N)(N^{-p/4}/\lambda)^{4/(p-2)}$, but the runtime increases
exponentially.

Restricting to the symmetric subspace, these ${n_{bos}}$ qudits can be thought
of as a system of bosons. In the case $p=4$, for example, our Hamiltonian has
pairwise interaction terms for all pairs of qudits. It is natural from the
viewpoint of mean-field theory in physics then to expect that the leading
eigenvector of the problem, for large ${n_{bos}}$, can be approximated by a
product state. While the bounds for arbitrary pairwise Hamiltonians would
require rather large ${n_{bos}}$ for given $N$ in order for such a mean-field
approximation to be accurate [8, 9, 10], we will be able to prove that for the
statistical model above the mean-field approximation becomes accurate with
high probablity at much smaller ${n_{bos}}$, depending upon the value of
$\lambda$. In this mean-field regime, the product state is an ${n_{bos}}$-fold
tensor product of a single particle state, and this single particle state is
given by $v_{\rm sig}$ in an obvious way, regarding the vector $v_{\rm sig}$
as a vector in the single-particle Hilbert space. While we will not prove that
this state is a good approximation to the leading eigenvector, it will be a
good approximation to some state in an eigenspace with large eigenvalue. Then,
the single particle density matrix allows one to infer $v_{\rm sig}$ (a
similar matrix was used in Ref. [1] where it was termed a voting matrix).

Classically, implementing this spectral algorithm requires high-dimensional
linear algebra, in particular finding the leading eigenvector of a matrix of
dimension $\approx N^{n_{bos}}$. This makes it a natural candidate for a
quantum algorithm. Since the Hamiltonian here is fairly simple, it can be
simulated efficiently using standard techniques in the literature reviewed
later. This allows us to give a simple algorithm based on preparing a random
initial state and then phase estimating in an attempt to project onto the
leading eigenvector. The probability of success in this projection is inverse
in the dimension of the matrix, so this simple algorithm leads to no speedup
over classical. However, we show that it is possible to apply amplitude
amplification to give a quadratic speedup over classical. More surprisingly,
we show that one can use the tensor $T_{0}$ to prepare an input state to the
algorithm with improved overlap with the leading eigenvector, giving the
quantum algorithm a quartic speedup over classical. Here, when comparing to
classical we are considering classical algorithms based on the power method or
similar algorithms such as Lanczos; these algorithms require exponential space
while the quantum algorithm uses only polynomial space. We also consider
classical algorithms based on ideas in Ref. [11] which use polynomial space
but the quantum algorithm is super-polynomially faster than these algorithms.
We also present some minor improvements to the quantum algorithm which may be
useful in practice.

###  1.1 Definitions, Random Ensembles, and Notation

Let us make some formal definitions. A tensor $T$ of order $p$ and dimension
$N$ is a multi-dimensional array. The entries of the tensor are written
$T_{\mu_{1},\mu_{2},\ldots,\mu_{p}}$ where $p\geq 1$ is an integer and each
$\mu_{a}$ ranges from $1,\ldots,N$. Generalizing previous work on this
problem, we consider two possible cases, one in which entries of a tensor are
chosen to be real numbers, and one in which they may be complex numbers, so
that either $T\in({\mathbb{R}^{N}})^{\otimes p}$ or
$T\in({\mathbb{C}}^{N})^{\otimes p}$; we explain later the reason for this
generalization; a tensor with all entries real will be called a real tensor. A
symmetric tensor is one that is invariant under permutation of its indices.
The symmetrization of a tensor is equal to $1/p!$ times the sum of tensors
given by permuting indices.

The spiked tensor model for given $N,p$ is defined as follows. Let $v_{\rm
sig}$ be a vector in ${\mathbb{R}}^{N}$, normalized by $|v_{\rm
sig}|=\sqrt{N}$, chosen from some probability distribution; this is the
“signal vector”. Let $G$ be a real tensor of order $p$ with entries chosen
from a Gaussian distribution with vanishing mean. We let $T_{0}=\lambda v_{\rm
sig}^{\otimes p}+G$ as above, where $v_{\rm sig}^{\otimes p}$ is defined to be
the tensor with entries

$(v_{\rm sig}^{\otimes p})_{\mu_{1},\ldots,\mu_{p}}=\prod_{a=1}^{p}(v_{\rm
sig})_{\mu_{a}}.$

Here we use the notation that a subscript on a vector denotes an entry of that
vector; we use a similar notation for matrices later.

Remark: some of the best sum-of-squares results are for a different
distribution in which the entries of $T_{0}$ are chosen from a biased
distribution on $\\{-1,+1\\}$, rather than for the Gaussian distribution. We
expect that using that distribution would not affect the results here too
much, but we avoid treating that case also for simplicity.

Since the tensor $v_{\rm sig}^{\otimes p}$ is symmetric, of course it is
natural to replace $T_{0}$ by its symmetrization. Indeed, no information can
be lost by this replacement since given a tensor $T_{0}$ one can symmetrize
the tensor, and then add back in Gaussian noise chosen to vanish under
symmetrization to obtain a tensor drawn from the same distribution as $T_{0}$
was. That is, the cases in which $G$ is symmetrized or not can be reduced to
each other.

A generalization of this problem is the case in which $G$ is chosen to have
complex entries, with each entry having real and imaginary parts chosen from a
Gaussian distribution with vanishing mean and variance $1/2$. We refer to this
as the complex ensemble, while we refer to the case where $G$ has real entries
as the real ensemble; the choice of reducing the variance to $1/2$ is a
convenient normalization for later. It is clear that since $v_{\rm sig}$ is
real, the case of complex $G$ can be reduced to the real case (up to an
overall rescaling for the different variance) simply by taking the real part
of $T_{0}$, and similarly the real case can be reduced to the complex case
(again up to an overall rescaling) by adding Gaussian distributed imaginary
terms to the entries of $T_{0}$. We will see that for odd $p$, at least for
reasons of analyzing the algorithms, it is convenient not to symmetrize
$T_{0}$ and to take complex $G$, while for even $p$ this is not necessary. It
may be possible to avoid doing this for odd $p$ (which may improve the
detection and recovery threshold of the algorithm by constant factors) and we
comment on this later.

We treat $p$ as fixed in the asymptotic notation, but consider the dependence
on $N,{n_{bos}}$. So, throughout this paper, when we refer to a polynomial in
$N$, we mean a polynomial independent of the parameter ${n_{bos}}$. The
polynomial may, however, depend on $p$, such as $N^{p}$.

We make additionally the following assumptions:

######  Assumption 1.

We assume that ${n_{bos}}=O(N^{\theta})$ for some $p$-dependent constant
$\theta>0$ chosen sufficiently small. We will also assume that $\lambda$ is
$\Omega(N^{-\theta^{\prime}})$ for some $p$-dependent constant
$\theta^{\prime}>p/4$.

Finally, we assume that $\lambda=O(N^{-p/4})$. Remark: there is of course no
reason to consider $\lambda$ larger than this since simple spectral methods
succeed if $\lambda$ is $\omega(N^{-p/4})$, but we state this assumption
explicitly as it simplifies some of the big-O notation.

We will explicitly state this Assumption 1 in all theorems where it is needed;
the assumption will be implicit in the statement of the lemmas and will not be
explicitly stated to avoid cluttering the statement of the results.

The first of these assumptions, that ${n_{bos}}=O(N^{\theta})$, is useful to
simplify some of the statements of the results to avoid having to specify the
allowed range of ${n_{bos}}$ in each case. For example, we will say that a
quantity such as ${n_{bos}}^{p}/N$ is $o(1)$, meaning that we must take
$\theta<1/p$. We do not specify the exact value of $\theta$ but it can be
deduced from the proofs if desired.

The second of these assumptions, that $\lambda$ is
$\Omega(N^{-\theta^{\prime}})$, also helps simplify some of the statements of
the results. Since we have assumed that ${n_{bos}}=O(N^{\theta})$ and we will
see that the required ${n_{bos}}$ increases polynomially with $\lambda^{-1}$,
this assumed lower bound on $\lambda$ is not a further restriction on our
results.

We write $\mathbb{E}[\ldots]$ to denote expectation values and ${\rm
Pr}[\ldots]$ to denote a probability. Usually these are expectation values or
probabilities over choices of $G$, though in some cases we consider
expectation values over other random variables. We use $\|\ldots\|$ to denote
the operator norm of an operator, i.e., the largest singular value of that
operator. We use $|\ldots|$ to denote either the $\ell_{2}$ norm of a tensor
or the $\ell_{2}$ norm of a quantum state. All logarithms are natural
logarithms.

We use $\langle\ldots|\ldots\rangle$ to denote inner products and use bra-ket
notation both for vectors and for quantum mechanical states.

We will need to compute expectation values over random choices of $G$ and also
compute expectation values of certain operators in quantum states, such as
$\langle\psi|O|\psi\rangle$ for some state $\psi$ and operator $O$. We refer
to the latter as a quantum expectation value of $O$ in state $\psi$ to
distinguish it from an expectation value over random variables.

###  1.2 Outline

In section 2, we review some results on recovery and boosting from Ref. [1]
and present a randomized recovery procedure that will help in improving the
recovery threshold. In section 3, we give spectral algorithms for the spiked
tensor problem for the case of both even and odd $p$. In that section, we
present algorithms in terms of eigenvalues and eigenvectors of a matrix (more
precisely, vectors in some eigenspace and quantum expectation values of
operators in those vectors) that we call a Hamiltonian. We leave the method to
compute these eigenvalues and expectation values for later in section 5, where
we give classical and quantum algorithms for this and give time bounds for
those algorithms. In section 4, we give some results on the spectra of random
tensors needed for the analysis of these algorithms. A key idea here is
reducing the case of a $p$-th order tensor for odd $p$ to the case of a $q$-th
order tensor for even $q=2(p-1)$. One interesting corollary of this technique,
see corollary 2, is that for odd $p$ and for the minimal value of ${n_{bos}}$
we are able remove a logarithmic factor in some of the bounds (a similar
logarithmic factor has been removed also using in [4] using what they termed
an “unfolding” algorithm). An appendix A gives an introduction to some
techniques used to evaluate expectation values of tensors networks whose
tensors are chosen from a Gaussian distribution; these techniques are used
earlier in the paper. In section 6, we further discuss tensor networks and use
this to consider limitations of certain algorithms and also to explain further
some of the motivation for this algorithm. In section 7, we discuss some
extensions of the results.

The proof of detection is in theorem 2 for the even $p$ case and 4 for the odd
$p$ case. The proof of recovery is in theorem 3 for the even $p$ case and
theorem 5 for the odd $p$ case. The runtime bound for the fastest quantum
algorithm is in theorem 6. This theorem gives a quartic improvement in the
runtime compared to the fastest classical spectral algorithm; more precisely
the log of the runtime with the quantum algorithm divided by the log of the
runtime of the classical algorithm approaches $1/4$ as $N\rightarrow\infty$ at
fixed $N^{-p/4}/\lambda$.

##  2 Recovery

In this section we discuss recovery and define randomized procedures for
recovery that will be useful in boosting the threshold. Following the notation
of Ref. [1], define the correlation between two vectors by a normalized
overlap

${\rm corr}(x,y)=\frac{|\langle x|y\rangle|}{|x|\cdot|y|}.$ (2)

The goal of an algorithm is to produce a vector $x$ with large ${\rm
corr}(x,v_{\rm sig})$. Note that we take the absolute value in the
correlation, ignoring the sign. For even $p$, the sign of $v_{\rm sig}$ is
irrelevant, while for odd $p$ it is easy, given a guess of $v_{\rm sig}$ up to
sign, to try both choices of sign and see which is most likely.

Strong recovery means that ${\rm corr}(x,v_{\rm sig})=1-o(1)$. Proposition 2.6
of Ref. [1], which is noted in that reference as being implicit in Ref. [2],
shows how to “boost” a weaker correlation to strong recovery. It is shown that
given a vector $u$ one can apply a single iteration of the the tensor power
algorithm to obtain a new vector $x$ such that, with high probability,

${\rm corr}(x,v_{\rm sig})\geq 1-c\lambda^{-1}{\rm corr}(u,v_{\rm
sig})^{1-p}N^{(1-p)/2},$ (3)

where $c$ is a constant depending on $p$. So, for any
$\lambda\omega(N^{(1-p)/2}$, if ${\rm corr}(u,v_{\rm sig})=\Omega(1)$, we have
${\rm corr}(x,v_{\rm sig})=1-o(1)$.

Thus, it suffices to construct which outputs some vector $u$ which has, with
high probability, ${\rm corr}(u,v_{\rm sig})=\Omega(1)$. This is termed weak
recovery; indeed, one can be satisfied with even weaker assumptions depending
on the value of $\lambda$; for $\lambda$ close to $N^{-p/4}$, one may have
polynomially small ${\rm corr}(u,v_{\rm sig})$.

The spectral algorithms that we construct later will output a matrix that we
will write $\rho_{1}$. This matrix will be a positive semi-definite Hermitian
matrix with trace $1$. In physics terms, such a matrix is termed a density
matrix. For sufficiently large $\lambda$, the leading eigenvector of the
matrix will have a large overlap with $v_{\rm sig}$. However, for smaller
$\lambda$, we will still have some lower bound that, with high probability,
$\langle v_{\rm sig}|\rho_{1}|v_{\rm sig}\rangle=\Omega(1)|v_{\rm sig}|^{2}$.
The randomized algorithm 1 below shows how to use this to obtain weak
recovery. This randomized algorithm allows us to improve the threshold over
the recovery method [1] of simply using the leading eigenvector $\rho_{1}$
since it works even in some cases where the leading eigenvector has small
correlation with $v_{\rm sig}$.

Putting these results together we find that

######  Corollary 1.

Given an algorithm that, with high probability, outputs a density matrix
$\rho$ with

$\frac{\langle v_{\rm sig}|\rho_{1}|v_{\rm sig}\rangle}{N}=\Omega(1),$

then in polynomial time, with high probability, strong recovery is possible.

We will simply say that such an algorithm achieves recovery.

We present the algorithm for the case that the matrix $\rho_{1}$ may be
complex; if the matrix is real, one can instead sample from a real Gaussian
distribution and the proof of lemma 1 goes through will slightly different
constants.

Algorithm 1 Input: density matrix $\rho$. Output, some vector $w$ obeying
bounds described above

  * 1.

Randomly sample a vector $u$ with entries chosen from a correlated complex
Gaussian distribution with zero mean and with covariance
$\mathbb{E}[\overline{u}_{i}u_{j}]=(\rho)_{ij},$ with
$\mathbb{E}[u_{i}u_{j}]=\mathbb{E}[\overline{u}_{i}\overline{u}_{j}]=0$.

  * 2.

Let $w=u/|u|$.

We have

######  Lemma 1.

For algorithm 1, with probability at least $1/2$,

$|\langle w|v_{\rm sig}\rangle|\geq c^{\prime}\sqrt{\langle v_{\rm
sig}|\rho|v_{\rm sig}\rangle},$ (4)

for some scalar $c^{\prime}>0$.

###### Proof.

We have $\mathbb{E}[|u|^{2}]={\rm tr}(\rho)=1$. Hence, with probability at
least $3/4$ we have that $|u|^{2}\leq 4/3$. We have $\mathbb{E}[|\langle
u|v_{\rm sig}\rangle|^{2}]=\langle v_{\rm sig}|\rho v_{\rm sig}\rangle$ and
since $\langle u|v_{\rm sig}\rangle$ is a Gaussian random variable with mean
$0$, with probability at least $3/4$ its absolute value is at least some
positive constant $c^{\prime\prime}$ (the exact constant can be deduced from
the error function) times its standard deviation. Hence, the lemma follows for
$c^{\prime}=(3/4)c^{\prime\prime}$. ∎

##  3 Spectral Algorithm

We now give the spectral algorithms for the spiked tensor problem for the case
of both even and odd $p$. In subsection 3.1, we define a Hamiltonian $H(T)$
given an arbitrary tensor $T$. Then, in subsection 3.2, we present the
spectral algorithm in terms of $H(T_{0})$.

For even $p$, the Hamiltonian that we present is very similar to the matrix
$Y$ given in Ref. [1] but it has some minor differences. In our language
(explained below), this matrix $Y$ is obtained by projecting our Hamiltonian
of Eq. (5) into the subspace of the “symmetric subspace” spanned by
$|\mu_{1}\rangle\otimes|\mu_{2}\rangle\otimes\ldots\otimes|\mu_{{n_{bos}}}\rangle$
with $\mu_{1},\ldots,\mu_{{n_{bos}}}$ all distinct. The relative reduction in
the size of the matrix is only $O(1/N)$ in the limit of large $N$.

Also, in our method, we have an $O(N)$ rotational symmetry of the basis which
is very useful in analysis, for example showing that the eigenvalues of
$H(\lambda v_{\rm sig}^{\otimes p})$ are independent of choice of $v_{\rm
sig}$. For the matrix $Y$ of [1], this is not obvious to us and we do not
fully understand the claimed behavior of the largest eigenvalue in that case.
We will use a different notation, using creation and annihilation operators,
which will help make this rotational symmetry more explicit.

For odd $p$, the Hamiltonian that we use is unrelated to that of Ref. [1].

###  3.1 Hamiltonian Definition

For even $p$, given a tensor $T$ we define an linear operator $H(T)$ that we
call a Hamiltonian as follows. This Hamiltonian is a linear operator on a
vector space $({\mathbb{R}}^{N})^{\otimes{n_{bos}}}$ or
$({\mathbb{C}}^{N})^{\otimes{n_{bos}}}$, for integer ${n_{bos}}\geq 1$ chosen
below. We write basis elements of this space as
$|\mu_{1}\rangle\otimes|\mu_{2}\rangle\otimes\ldots\otimes|\mu_{n_{bos}}\rangle$,
and we call this space the full Hilbert space. We define

$H(T)=\frac{1}{2}\sum_{i_{1},\ldots,i_{p/2}}\Bigl{(}\sum_{\mu_{1},\ldots,\mu_{p}}T_{\mu_{1},\mu_{2},\ldots,\mu_{p}}|\mu_{1}\rangle_{i_{1}}\langle\mu_{1+p/2}|\otimes|\mu_{2}\rangle_{i_{2}}\langle\mu_{2+p/2}|\otimes\ldots\otimes|\mu_{p/2}\rangle_{i_{p/2}}\langle\mu_{p}|+{\rm
h.c.}\Bigr{)},$ (5)

where the sum is over distinct $i_{1},i_{2},\ldots,i_{p/2}$ so that there are
$(p/2!){{n_{bos}}\choose p/2}$ terms in the sum and where ${\rm h.c.}$ means
adding the Hermitian conjugate of the given terms, so that $H(T)$ is Hermitian
and where $|\mu\rangle_{i}\langle\nu|$ denotes the operator
$|\mu\rangle\langle\nu|$ on qudit $i$. We require that ${n_{bos}}\geq p/2$ or
else $H(T)$ is trivially zero.

Note of course that if $T$ is real and symmetric, then the term
$\sum_{\mu_{1},\ldots,\mu_{p}}T_{\mu_{1},\mu_{2},\ldots,\mu_{p}}|\mu_{1}\rangle_{i_{1}}\langle\mu_{1+p/2}|\otimes|\mu_{2}\rangle_{i_{2}}\langle\mu_{2+p/2}|\otimes\ldots\otimes|\mu_{p/2}\rangle_{i_{p/2}}\langle\mu_{p}|$
is already Hermitian. $H(T)$ can be regarded as a Hamiltonian acting on a
space of ${n_{bos}}$ qudits, each of dimension $N$, and with interaction
between sets of $p/2$ particles at a time.

Even if $T$ is not symmetrized, $H(T)$ is unchanged if one applies an
arbitrary permutation to the first $p/2$ indices of $T$ and applies the same
permutation to the last $p/2$ indices of $T$.

We may restrict to the symmetric subspace of this Hilbert space. We write
$D(N,{n_{bos}})$ to indicate the dimension of this subspace. For
$N\gg{n_{bos}}$, we can approximate $D(N,{n_{bos}})\approx
N^{n_{bos}}/{n_{bos}}!$.

Within the symmetric subspace, we can write this Hamiltonian in a so-called
“second-quantized” form:

$H(T)=\frac{1}{2}\Bigl{(}\sum_{\mu_{1},\ldots,\mu_{p}}T_{\mu_{1},\mu_{2},\ldots,\mu_{p}}\Bigl{(}\prod_{i=1}^{p/2}a^{\dagger}_{\mu_{i}}\Bigr{)}\Bigl{(}\prod_{i=p/2+1}^{p}a_{\mu_{i}}\Bigr{)}+{\rm
h.c.}\Bigr{)}.$ (6)

This replacement by a second-quantized Hamiltonian is simply a convenient
notation. The operators $a^{\dagger}_{\mu},a_{\mu}$ are bosonic creation and
annihilation operators, obeying canonical commutation relations
$[a_{\mu},a^{\dagger}_{\nu}]=\delta_{\mu,\nu}$. See appendix B for a brief
review of this formalism. We restrict to the subspace with a total of
${n_{bos}}$ bosons, i.e., we define the number operator $n$ by

$n\equiv\sum_{\mu}a^{\dagger}_{\mu}a_{\mu},$ (7)

and restrict to $n={n_{bos}}.$ An orthonormal basis of states for this
symmetric subspace is given by all states equal to some normalization constant
multiplying $a^{\dagger}_{\mu_{1}}a^{\dagger}_{\mu_{2}}\ldots
a^{\dagger}_{\mu_{n_{bos}}}|0\rangle$, where $|0\rangle$ is the vacuum state
(i.e., the state annihilated by $a_{\mu}$ for all $\mu$), and where
$\mu_{1}\leq\mu_{2}\leq\ldots\leq\mu_{n_{bos}}$ is some sequence.

The second quantized Hamiltonian for the symmetric subspace is unchanged under
arbitrary permutation of the first $p/2$ indices of $T$ and arbitrary (not
necessarily the same) permutation of the last $p/2$ indices of $T$.

For odd $p$, we define the Hamiltonian $H(T)$ as follows. Given a tensor $T$
of odd order $p$, define a new tensor $\tilde{T}$ of even order $q=2(p-1)$
with components

$\tilde{T}_{\mu_{1},\ldots,\mu_{(p-1)/2},\nu_{1},\ldots,\nu_{(p-1)/2},\mu_{(p-1)/2+1},\ldots,\mu_{p-1},\nu_{(p-1)/2},\ldots,\nu_{p-1}}=\sum_{\sigma}T_{\mu_{1},\ldots,\mu_{p-1},\sigma}T_{\nu_{1},\ldots,\nu_{p-1},\sigma}.$
(8)

Then define $H(T)=H(\tilde{T})$, using the definition (5) for $H(\tilde{T})$.
Note the order of indices on the left-hand side of Eq. (8). Using the second-
quantized notation, this gives for odd $p$:

$\displaystyle H(T)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\Bigl{(}\sum_{\mu_{1},\ldots,\mu_{p-1}}\sum_{\nu_{1},\ldots,\nu_{p-1}}\sum_{\sigma}T_{\mu_{1},\mu_{2},\ldots,\mu_{p-1},\sigma}T_{\nu_{1},\nu_{2},\ldots,\nu_{p-1},\sigma}\Bigl{(}\prod_{i=1}^{(p-1)/2}a^{\dagger}_{\mu_{i}}a^{\dagger}_{\nu_{i}}\Bigr{)}\Bigl{(}\prod_{i=(p-1)/2+1}^{p-1}a_{\mu_{i}}a_{\nu_{i}}\Bigr{)}+{\rm
h.c.}\Bigr{)},$

Now we require that ${n_{bos}}\geq p-1$ as otherwise $H(T)$ is trivially zero.
For this Hamiltonian, it is convenient to take $G$ from the complex ensemble
because, as we explain more below, it makes $\mathbb{E}[H(G)]$ equal to zero,
as well as canceling out certain terms in higher order moments, making the
proof of the spectral properties of $H(G)$ simpler. We discus later to what
extent we can avoid using the complex ensemble.

###  3.2 Spectral Algorithms

The spectral algorithm for detection and recovery is algorithm 2. In this
subsection we prove correctness of this algorithm, using statistical
properties of $H(G)$ proven later.

This algorithm uses quantities $E_{0}$ and $E_{max}$ defined later; roughly
$E_{max}$ is an upper bound on the eigenvalues of $H(G)$ and $E_{0}$ is the
largest eigenvalue of $H(\lambda v^{\otimes p})$. The algorithm can then
achieve detection by verifying that the largest eigenvalue is significantly
larger than $E_{max}$, which occurs when $E_{0}$ is large enough. Indeed, we
will see that it suffices to have $E_{0}=(1+c)E_{max}$ for any fixed $c>0$
(some results can be extended to the case that $c$ decays slowly with $N$ but
we omit this for brevity). This fixes the scaling of ${n_{bos}}$ with
$\lambda$ so that we need (up to polylogarithmic factors)
${n_{bos}}\gtrsim(N^{-p/4}/\lambda)^{4/(p-2)}$.

One interesting feature of the algorithm is that in step 3, we compute the
density matrix of the leading eigenvector or of any vector in the eigenspace
of eigenvalue $\geq E_{cut}$, for $E_{cut}=(E_{0}+E_{max})/2$ defined in the
algorithm. This might seem surprising, given that the leading eigenvector is
computed in step 1; one might wonder why some other vector should be taken. We
describe the algorithm in this way since, in later classical and quantum
algorithms that we give to compute the spectral properties of the matrix, we
might not extract the leading eigenvector but instead extract only some vector
in this eigenspace due to use of the power method in a classical algorithm or
due to approximations in phase estimation in a quantum algorithm. Thus, much
of our analysis is given to showing that some eigenvalue larger than $E_{cut}$
exists by lower bounding the leading eigenvalue $\lambda_{1}$, but given that
some such eigenvalue exists, we do not worry too much about exactly what
mixture of eigenvectors in the given eigenspace we compute.

Algorithm 2 Spectral algorithm. This algorithm takes a tensor $T_{0}$ as input
and also a scalar $\overline{\lambda}$ and an integer ${n_{bos}}$. The output
is a decision about whether $\lambda=\overline{\lambda}$ or $\lambda=0$, and,
if the algorithm reports that $\lambda=\overline{\lambda}$, it also returns an
approximation of $v_{\rm sig}$ (up to an overall sign). The quantity
${n_{bos}}$ is chosen depending upon the value of $\overline{\lambda}$;
smaller values of $\lambda$ require larger values of ${n_{bos}}$ in order for
$E_{0}$ to be sufficiently larger than $E_{max}$ for the algorithm to be
accurate. See theorems 2,4,3,5. For $E_{0}\geq(1+c)E_{max}$ for any $c>0$, the
algorithm achieves recovery.

  * 1.

Compute the eigenvector of $H(T_{0})$ and the leading eigenvalue, denoted
$\lambda_{1}$.

  * 2.

(Detection) If

$\lambda_{1}>E_{cut}\equiv(E_{0}+E_{max})/2,$

where $E_{0}=\overline{\lambda}(p/2)!{{n_{bos}}\choose p/2}N^{p/2}$ for even
$p$, and $E_{0}=\overline{\lambda}^{2}(p-1)!{{n_{bos}}\choose p-1}N^{p}$ for
odd $p$, and where $E_{max}$ is defined in theorem 1, then report that
$\lambda=\overline{\lambda}$. Otherwise report $\lambda=0$.

  * 3.

(Recovery) Compute the single particle density matrix (defined below) of the
leading eigenvector or of any vector in the eigenspace of eigenvalue $\geq
E_{cut}$. Apply algorithm 1 to recover an approximation to $v_{\rm sig}$.

In section 4, we consider the largest eigenvalue $\lambda_{1}$ of $H(G)$ and
show the following theorem which summarizes the results of lemmas 5,7.
Roughly, up to prefactors, the result for odd $p$ is given by considering the
result for a $q$-th order tensor for even $q=2(p-1)$ and then multiplying the
largest eigenvalue by a factor of $\sqrt{N}$.

######  Theorem 1.

Let $\lambda_{1}$ be the largest eigenvalue of $G$. For even $p$ let

$E_{max}=\sqrt{2J\log(N)}{n_{bos}}^{p/4+1/2}N^{p/4},$

and for odd $p$ let

$E_{max}=2\sqrt{J\log(N)}{n_{bos}}^{p/2}N^{p/2},$

where $J$ is a scalar depends that implicitly on $p,{n_{bos}},N$ and tends to
some function depending only on $p$ for large ${n_{bos}},N$. More precisely,
for even $p$, $J$ is equal to $(p/2)!{{n_{bos}}\choose
p/2!}/{n_{bos}}^{p/2}+o(1)$ for the real ensemble and is twice that for the
complex ensemble, and for odd $p$, $J$ is equal to that for the even case for
$2(p-1)$.

Then, for any $x$, assuming Assumption 1,

${\rm Pr}[\lambda_{1}\geq x]\leq\exp\Bigl{(}-\frac{x-E_{max}}{\xi}\Bigr{)},$
(10)

with for even $p$

$\xi=\frac{\sqrt{J}{n_{bos}}^{p/4-1/2}N^{p/4}}{\sqrt{2\log(N)}}$ (11)

and for odd $p$

$\xi=\frac{\sqrt{J}{n_{bos}}^{p/2-1}N^{p/2}}{\sqrt{\log(N)}}.$ (12)

So, for any $E^{\prime}$ which is $\omega(\xi)$, with high probability
$\lambda_{1}\leq E_{max}+E^{\prime}$.

Now consider the eigenvectors and eigenvalues of $H(T_{0})$. For any symmetric
tensor $T$ of order ${n_{bos}}$, let $|T\rangle$ be the vector on ${n_{bos}}$
qudits (each of dimension $N$) with amplitudes given by the entries of the
tensor in the obvious way:

$|T\rangle=\sum_{\mu_{1},\ldots,\mu_{{n_{bos}}}}T_{\mu_{1},\ldots,\mu_{{n_{bos}}}}|\mu_{1}\rangle\otimes\ldots\otimes|\mu_{{n_{bos}}}\rangle.$

This vector is only normalized if $|T|=1$. So, $\Psi_{\rm sig}\equiv
N^{-{n_{bos}}/2}|v_{\rm sig}^{\otimes{n_{bos}}}\rangle$ is a normalized
vector.

We have the following simple property:

######  Lemma 2.

Let $\lambda_{1}$ be the largest eigenvalues of $H(T_{0})$. Then,
$\lambda_{1}\geq E\equiv\langle\Psi_{\rm sig}|H(T)|\Psi_{\rm sig}\rangle$.

###### Proof.

Immediate from the variational principle. ∎

####  3.2.1 Even $p$ Case

We now show correctness of the algorithm. All results in this subsubsection
refer to even $p$ even if we do not state it explicitly. First we estimate $E$
of lemma 2 to show detection. Then, we show recovery.

We have

$\displaystyle E$ $\displaystyle=$ $\displaystyle\langle\Psi_{\rm
sig}|H(\lambda v_{\rm sig}^{\otimes p})|\Psi_{\rm sig}\rangle+\langle\Psi_{\rm
sig}|H(G)|\Psi_{\rm sig}\rangle$ (13) $\displaystyle=$ $\displaystyle
E_{0}+\langle\Psi_{\rm sig}|H(G)|\Psi_{\rm sig}\rangle,$ (14)

where

$E_{0}=\lambda(p/2)!{{n_{bos}}\choose p/2}N^{p/2}.$ (15)

To evaluate $\langle\Psi_{\rm sig}|H(G)|\Psi_{\rm sig}\rangle$ it is
convenient to exploit a rotational invariance of this problem. We can apply a
rotation using any matrix $O$ in the orthogonal group $O(N)$, rotating the
creation and annihilation operators by making the replacement

$a^{\dagger}_{\mu}\rightarrow\sum_{\nu}O_{\mu\nu}a^{\dagger}_{\nu}$

and

$a_{\mu}\rightarrow\sum_{\nu}O_{\mu\nu}a_{\nu}.$

This rotation preserves the canonical commutation relations and is equivalent
to rotating the basis states on each qudit by $O$. To preserve the Hamiltonian
$H(T_{0})$, we rotate each leg of the tensor $T_{0}$:

$(T_{0})_{\mu_{1},\ldots,\mu_{p}}\rightarrow\sum_{\nu_{1},\ldots,\nu_{p}}(T_{0})_{\nu_{1},\ldots,\nu_{p}}\prod_{i}O_{\mu_{i},\nu_{i}}.$

This rotation preserves the Gaussian measure on $G$ but changes $v_{\rm sig}$.
So, we can rotate so that $v_{\rm sig}$ is some fixed vector, say
$(\sqrt{N},0,0,\ldots)$ so $(v_{\rm sig})_{1}=\sqrt{N}$. Then
$\langle\Psi_{\rm sig}|H(G)|\Psi_{\rm sig}\rangle$ is equal to
$(p/2)!{{n_{bos}}\choose p/2}$ multiplied by a single entry of $G$, i.e., the
entries with all indices equal to $1$, which is some quantity chosen from a
Gaussian distribution with zero mean and unit variance. So, with high
probability, $E=\lambda(p/2)!{{n_{bos}}\choose
p/2}N^{p/2}+{n_{bos}}^{p/2}O(N^{\eta})$ for any $\eta>0$.

Hence,

######  Theorem 2.

If $\lambda(p/2)!{{n_{bos}}\choose
p/2}N^{p/2}-E_{max}=\omega\Bigl{(}\frac{\sqrt{J}{n_{bos}}^{p/4-1/2}N^{p/4}}{\sqrt{2\log(N)}}\Bigr{)}$,
and if Assumption 1 holds, then with high probability algorithm 2 correctly
determines whether $\lambda=0$ or $\lambda=\overline{\lambda}$.

###### Proof.

This follows from the estimate of $E$ which lowers bounds the largest
eigenvalue when $\lambda=\overline{\lambda}$ and from theorem 1 which upper
bounds $\|H(G)\|$. ∎

Given any state, define the single particle density matrix $\rho_{1}$ in the
basis for the full Hilbert space used in Eq. (5) to be the reduced density
matrix on any of the qudits. Equivalently, this single particle density matrix
can be expressed in terms of creation and annihilation operators by

$(\rho_{1})_{\mu\nu}=\frac{1}{{n_{bos}}}a^{\dagger}_{\mu}a_{\nu}.$ (16)

Note that $\rho_{1}$ is positive semi-definite, Hermitian, and trace $1$.

We have

######  Theorem 3.

Let Assumption 1 hold. Given any vector $\Psi$ such that
$\langle\Psi|H(T_{0})|\Psi\rangle\geq(1+c^{\prime})E_{max}$ for any scalar
$c^{\prime}>0$, then with high probability the corresponding single particle
density matrix $\rho_{1}$ obeys

$\frac{\langle v_{\rm sig}|\rho_{1}|v_{\rm
sig}\rangle}{N}\geq(c^{\prime}-o(1))\frac{E_{max}}{E_{0}}.$ (17)

In particular, for $\langle\Psi|H(T_{0})|\Psi\rangle\geq E_{cut}$ with
$E_{0}\geq(1+c)E_{max}$ we have $c^{\prime}\geq c/2$ so with high probability

$\frac{\langle v_{\rm sig}|\rho_{1}|v_{\rm
sig}\rangle}{N}\geq\frac{1}{2}\frac{c}{1+c}.$ (18)

Hence, algorithm 2 achieves recovery.

###### Proof.

We have $\langle\Psi|H(\lambda v_{\rm sig}^{\otimes
p}|\Psi\rangle+\langle\Psi|H(G)|\Psi\rangle\geq(1+c)E_{max}$. By theorem 1,
with high probability $\langle\Psi|H(G)|\Psi\rangle\leq(1+o(1))E_{max}$.
Hence, with high probability $\langle\Psi|H(\lambda v_{\rm sig}^{\otimes
p})|\Psi\rangle\geq(c-o(1))E_{max}$.

Rotate so that $v_{\rm sig}=(\sqrt{N},0,0,\ldots)$. Then, the Hamiltonian
$H(\lambda v_{\rm sig}^{\otimes p})$ is diagonal in the computational basis
for the full Hilbert space. Let $P(m)$ be the probability, for state $\Psi$,
that exactly $m$ out of ${n_{bos}}$ qudits are in state $|1\rangle$. Then,
$\langle\Psi|H(\lambda v_{\rm sig}^{\otimes
p})|\Psi\rangle=\sum_{m}\lambda(p/2)!{m\choose
p/2}N^{p/2}P(m)\leq\sum_{m}m(E_{0}/{n_{bos}})P(m)=E_{0}\langle v_{\rm
sig}|\rho_{1}|v_{\rm sig}\rangle/N$. ∎

Remark: more generally, the same result holds for any mixed state: given any
mixed state $\sigma$ such that ${\rm tr}(\sigma H(T_{0}))\geq(1+c)E_{max}$ for
any scalar $c>0$, then with high probability the corresponding single particle
density matrix $\rho_{1}$ obeys Eq. (17).

####  3.2.2 Odd $p$ Case

We now show correctness of the algorithm for odd $p$. All results in this
subsubsection refer to odd $p$ even if we do not state it explicitly.

Let us first estimate $E$ to show detection. We have $E=\langle\Psi_{\rm
sig}|H(T_{0})|\Psi_{\rm sig}\rangle$, but now $H(T_{0})$ is a quadratic
function of $T_{0}$ so there are some cross-terms. We have

$\langle\Psi_{\rm sig}|H(v_{\rm sig}^{\otimes p})|\Psi_{\rm
sig}\rangle=\lambda^{2}(p-1)!{{n_{bos}}\choose p-1}N^{p}\equiv E_{0}.$ (19)

The cross-term is

$\lambda\langle\Psi_{\rm
sig}|\sum_{\mu_{1},\ldots,\mu_{p-1}}\sum_{\nu_{1},\ldots,\nu_{p-1}}\sum_{\sigma}(v_{\rm
sig}^{\otimes}p)_{\mu_{1},\mu_{2},\ldots,\mu_{p-1},\sigma}G_{\nu_{1},\nu_{2},\ldots,\nu_{p-1},\sigma}\Bigl{(}\prod_{i=1}^{(p-1)/2}a^{\dagger}_{\mu_{i}}a^{\dagger}_{\nu_{i}}\Bigr{)}\Bigl{(}\prod_{i=(p-1)/2+1}^{p-1}a_{\mu_{i}}a_{\nu_{i}}\Bigr{)}+{\rm
h.c.}|\Psi_{\rm sig}\rangle.$

Exploiting the same rotational invariance as in the even $p$ case, this is
equal to $2\lambda(p-1)!{{n_{bos}}\choose p-1}N^{p/2}$ multiplied by the real
part of a single entry of $G$, i.e., the entry with all indices equal to $1$.
So, with high probability, this cross-term is bounded by
$\lambda{n_{bos}}^{p-1}O(N^{p/2+\eta})$ for any $\eta>0$.

The term quadratic in $G$ is

$\frac{1}{2}\langle\Psi_{\rm
sig}|\sum_{\mu_{1},\ldots,\mu_{p-1}}\sum_{\nu_{1},\ldots,\nu_{p-1}}\sum_{\sigma}G_{\mu_{1},\mu_{2},\ldots,\mu_{p-1},\sigma}G_{\nu_{1},\nu_{2},\ldots,\nu_{p-1},\sigma}\Bigl{(}\prod_{i=1}^{(p-1)/2}a^{\dagger}_{\mu_{i}}a^{\dagger}_{\nu_{i}}\Bigr{)}\Bigl{(}\prod_{i=(p-1)/2+1}^{p-1}a_{\mu_{i}}a_{\nu_{i}}\Bigr{)}+{\rm
h.c.}|\Psi_{\rm sig}\rangle.$

Exploiting the same rotational invariance as before, fixing $v_{\rm
sig}=(\sqrt{N},0,0,\ldots)$, we must have all $\mu_{i},\nu_{i}=1$. So, this is
a sum of squares of $N$ different entries of $G$, corresponding to $N$
possible choices of $\sigma$; since we use the complex ensemble, this has
vanishing mean. So, with high probability this term is
${n_{bos}}^{p-1}O(N^{1/2+\eta})$ for any $\eta>0$.

So, with high probability,

$E\geq
E_{0}-\lambda{n_{bos}}^{p-1}O(N^{p/2+\eta})-{n_{bos}}^{p-1}O(N^{1/2+\eta})$
(20)

for any $\eta>0$.

Hence,

######  Theorem 4.

Let Assumption 1 hold. If $\lambda^{2}(p-1)!{{n_{bos}}\choose
p-1}N^{p}-E_{max}=\omega\Bigl{(}\frac{\sqrt{J}{n_{bos}}^{p/2-1}N^{p/2}}{\sqrt{\log(N)}}\Bigr{)}$,
then with high probability algorithm 2 correctly determines whether
$\lambda=0$ or $\lambda=\overline{\lambda}$.

###### Proof.

This follows from the estimate of $E$ in Eq. (20) which lowers bounds the
largest eigenvalue when $\lambda=\overline{\lambda}$ and from theorem 1 which
upper bounds $\|H(G)\|$. Note that for $\lambda=O(N^{-p/4})$, the terms
$-\lambda{n_{bos}}^{p-1}O(N^{p/2+\eta})-{n_{bos}}^{p-1}O(N^{1/2+\eta})$ on the
right-hand side of this equation are negligible as they are
$o(\Bigl{(}\frac{\sqrt{J}{n_{bos}}^{p/2-1}N^{p/2}}{\sqrt{\log(N)}}\Bigr{)}$. ∎

We now consider recovery. We will first give a more general bound on cross-
terms that will be useful. We have

$H(T_{0})=H(\lambda v_{\rm sig}^{\otimes p})+H(G)+{\rm cross\,terms}.$

Let us first bound the cross-terms.

######  Lemma 3.

With high probability, the operator norm of the cross-terms is bounded by
$\lambda N^{p/2}O(N^{(p-1)/4}){n_{bos}}^{p-1}$.

###### Proof.

The cross-terms are equal to $2H(T_{cross})$ where $T_{cross}$ is a tensor for
even $q=2(p-1)$ which has components

$\Bigl{(}T_{cross}\Bigr{)}_{\mu_{1},\ldots,\mu_{(p-1)/2},\nu_{1},\ldots,\nu_{(p-1)/2},\mu_{(p-1)/2+1},\ldots,\mu_{p},\nu_{(p-1)/2},\ldots,\nu_{p}}=\lambda\sum_{\sigma}T_{\mu_{1},\ldots,\mu_{p-1},\sigma}\cdot\Bigl{(}\prod_{i=1}^{p-1}(v_{\rm
sig})_{\nu_{i}}\Bigr{)}\cdot(v_{\rm sig})_{\sigma}.$ (21)

Rotating so that $v_{\rm sig}=(\sqrt{N},0,0,\ldots)$, we have

$\Bigl{(}T_{cross}\Bigr{)}_{\mu_{1},\ldots,\mu_{(p-1)/2},\nu_{1},\ldots,\nu_{(p-1)/2},\mu_{(p-1)/2+1},\ldots,\mu_{p},\nu_{(p-1)/2},\ldots,\nu_{p}}=\lambda
N^{p/2}\sum_{\sigma}T_{\mu_{1},\ldots,\mu_{p-1},1}\prod_{i=1}^{p-1}\delta_{\nu_{i},1}.$
(22)

Clearly, $\|H(T_{cross})\|\leq{n_{bos}}^{p-1}\|M_{cross}\|,$ where $M_{cross}$
is an $N^{p-1}$-by-$N^{p-1}$ matrix, whose entries are determined by the
entries of $T_{cross}$ in an obvious way so that the first $p-1$ indices of
$T_{cross}$ determine a column index in $M$ and the last $p-1$ indices of
$T_{cross}$ determine a row index of $M$. Regard
$T_{\mu_{1},\ldots,\mu_{p-1},1}$ as being the entries of some tensor of order
$p-1$ and let $M^{\prime}$ by the $N^{(p-1)/2}$-by-$N^{(p-1)/2}$ matrix whose
entries are determined by the entries of this tensor, again in the obvious way
so that the first $(p-1)/2$ indices of the tensor determine a column index and
the last $(p-1)/2$ indices determine a row index.

Then,

$\|M_{cross}\|=\lambda N^{p/2}\|M^{\prime}\|.$ (23)

However, since the entries of $M^{\prime}$ are independent Gaussian random
variables, it follows from standard random matrix theory results [12] that
with high probability $\|M^{\prime}\|=O(N^{(p-1)/4})$. Remark: the dependence
on ${n_{bos}}$ is not tight here. ∎

So as in the even case, we have

######  Theorem 5.

Let Assumption 1 hold. Given any vector $\Psi$ such that
$\langle\Psi|H(T_{0})|\Psi\rangle\geq(1+c^{\prime})E_{max}$ for any scalar
$c^{\prime}>0$, then with high probability the corresponding single particle
density matrix $\rho_{1}$ obeys

$\frac{\langle v_{\rm sig}|\rho_{1}|v_{\rm
sig}\rangle}{N}\geq(c^{\prime}-o(1))\frac{E_{max}}{E_{0}}.$ (24)

In particular, for $\langle\Psi|H(T_{0})|\Psi\rangle\geq E_{cut}$ with
$E_{0}\geq(1+c)E_{max}$ we have $c^{\prime}\geq c/2$ so with high probability

$\frac{\langle v_{\rm sig}|\rho_{1}|v_{\rm
sig}\rangle}{N}\geq\frac{1}{2}\frac{c}{1+c}.$ (25)

Hence, algorithm 2 achieves recovery.

###### Proof.

We use lemma 3 to bound $\|H(T_{cross})\|$. This bound is asymptotically
negligible compared to $E_{max}$. So, we have $\langle\Psi|H(\lambda v_{\rm
sig}^{\otimes
p}|\Psi\rangle+\langle\Psi|H(G)|\Psi\rangle\geq(1+c-o(1))E_{max}$ where the
$o(1)$ denotes the $o(1)$ contribution from the cross-terms.

Then, the rest of the proof is the same as in the even case, except for a
replacement of $p/2$ by $p-1$. In detail: by theorem 1, with high probability
$\langle\Psi|H(G)|\Psi\rangle\leq(1+o(1))E_{max}$. Hence,
$\langle\Psi|H(\lambda v_{\rm sig}^{\otimes
p}|\Psi\rangle\geq(c-o(1))E_{max}$.

Rotate so that $v_{\rm sig}=(\sqrt{N},0,0,\ldots)$. Then, the Hamiltonian
$H(\lambda v_{\rm sig}^{\otimes p})$ is diagonal in the computational basis
for the full Hilbert space. Let $P(m)$ be the probability, for state $\Psi$,
that exactly $m$ out of ${n_{bos}}$ qudits are in state $|1\rangle$. Then,
$\langle\Psi|H(\lambda v_{\rm sig}^{\otimes
p})|\Psi\rangle=\sum_{m}\lambda^{2}(p-1)!{m\choose
p-1}N^{p}P(m)\leq\sum_{m}m(E_{0}/{n_{bos}})P(m)=E_{0}\langle v_{\rm
sig}|\rho_{1}|v_{\rm sig}\rangle/N$. ∎

##  4 Spectrum of Random Hamiltonian

In this section, we will estimate the eigenvalues of $H(G)$. We consider first
the case of even $p$. Here our proof is very similar to the of Ref. [1],
though the method here also suggests some heuristics that may lead to a
tighter bound in the future. Then, we consider the case of odd $p$ by reducing
it to the case of even $p$.

###  4.1 Even $p$

We first consider the case of even $p$. Let $Z(\tau,p,N,{n_{bos}})$ denote
$\mathbb{E}[{\rm tr}(\exp\\{\tau H(G)\\})]$ for a tensor $G$ of rank $p$, with
entries of $G$ chosen from the Gaussian distribution, with given
$N,{n_{bos}}$, for some real scalar $\tau$. In this subsection, $G$ may be
symmetrized or not, and may be chosen from either the real or complex
ensemble.

The main result in this section is the following lemma:

######  Lemma 4.

For each ensemble, real or complex, symmetrized or not, we have

$Z(\tau,p,N,{n_{bos}})\leq
D(N,{n_{bos}})\exp(\tau^{2}(J/2){n_{bos}}^{p/2}N^{p/2}),$ (26)

where $J$ is a scalar depends that implicitly on $p,{n_{bos}},N$ and tends to
some function depending only on $p$ for large ${n_{bos}},N$. More precisely,
$J$ is equal to $(p/2)!{{n_{bos}}\choose p/2!}/{n_{bos}}^{p/2}+o(1)$ for the
real ensemble and is twice that for the complex ensemble.

###### Proof.

We first give a brief derivation of Eq. (4) below which is a standard result
using quantum field theory techniques. Note that $\tau H(G)=H(\tau G)$ and
$\tau G$ is a tensor with entries chosen from a Gaussian distribution with
zero mean and variance $\tau^{2}$. Hence for any $\tau^{\prime}>\tau$ we have

$Z(\tau^{\prime},p,N,{n_{bos}})=\mathbb{E}_{G,\eta}[{\rm tr}(\exp\\{H(\tau
G+\eta)\\})],$ (27)

where $\mathbb{E}_{G,\eta}[\ldots]$ denotes the expectation value over $G$ and
$\eta$, with the tensor $\eta$ having Gaussian entries with zero mean and
variance $(\tau^{\prime})^{2}-\tau^{2}$. Taking the expectation value over
$\eta$, for $\tau^{\prime 2}-\tau^{2}$ small, we need to keep only the zeroth
and second order terms on the right-hand side. So, we find that

$\displaystyle\partial_{\tau}^{2}Z(\tau,p,N,{n_{bos}})$ $\displaystyle=$
$\displaystyle\int_{0}^{1}{\rm d}s_{2}\int_{0}^{s_{2}}{\rm
d}s_{1}\,\mathbb{E}_{G,\eta}[{\rm tr}\Bigl{(}\exp\Bigl{\\{}(1-s_{2})\tau
H(G)\Bigr{\\}}H(\eta)\exp\Bigl{\\{}(s_{2}-s_{1})\tau
H(G)\Bigr{\\}}H(\eta)\exp\Bigl{\\{}s_{1}\tau H(G)\Bigr{\\}}\Bigr{)}].$

Using cyclic properties of the trace, this can be simplified to

$\partial_{\tau}^{2}Z(\tau,p,N,{n_{bos}})=\frac{1}{2}\int_{0}^{1}{\rm
d}s_{1}\,\mathbb{E}_{G,\eta}[{\rm tr}\Bigl{(}\exp\Bigl{\\{}(1-s_{1})\tau
H(G)\Bigr{\\}}H(\eta)\exp\Bigl{\\{}s_{1}\tau H(G)\Bigr{\\}}H(\eta)\Bigr{)}].$
(29)

We now use a general result. Consider any Hermitian $H$ (we will use $H=\tau
H(G)$) and any operator $O$ (we will use $O=H(\eta)=H(\eta)^{\dagger}$) and
any $s_{1}\in[0,1]$. We claim that

${\rm tr}\Bigl{(}\exp(H)O^{\dagger}O\Bigr{)}+O\leftrightarrow
O^{\dagger}\geq{\rm
tr}\Bigl{(}\exp((1-s_{1})H)O^{\dagger}\exp(s_{1}H)O\Bigr{)}+O\leftrightarrow
O^{\dagger},$

where $O\leftrightarrow O^{\dagger}$ indicates the previous term with
$O,O^{\dagger}$ interchanged. Proof of claim: work in an eigenbasis of $H$. It
suffices to consider the case that $O=|b\rangle\langle a|$ where
$|a\rangle,|b\rangle$ are eigenvectors of $H$ with eigenvalues $E_{a},E_{b}$.
Then the the right-hand side is equal to
$\exp(s_{1}E_{a}+(1-s_{1})E_{b})+\exp(s_{1}E_{b}+(1-s_{1})E_{a})=2\exp((E_{a}+E_{b})/2)\cosh((s_{1}-1/2)(E_{b}-E_{a})$.
The cosh function is maximized on the interval $[0,1]$ at $s_{1}=0,1$ when the
right-hand side becomes equal to the left-hand side. So,

$\partial_{\tau}^{2}Z(\tau,p,N,{n_{bos}})\leq\frac{1}{2}\mathbb{E}_{G,\eta}[{\rm
tr}\Bigl{(}\exp(\tau H(G))H(\eta)^{2}\Bigr{)}].$ (30)

For the real ensemble without symmetrization, we have

$\mathbb{E}[H(\eta)^{2}]=\sum_{\mu_{1},\ldots,\mu_{p/2}}\sum_{\nu_{1},\ldots,\nu_{p/2}}\prod_{i=1}^{p/2}\Bigl{(}a^{\dagger}_{\mu_{i}}a_{\nu_{i}}\Bigr{)}\prod_{i=1}^{p/2}\Bigl{(}a^{\dagger}_{\nu_{i}}a_{\mu_{i}}\Bigr{)}\Bigr{.}$

To leading order in $N$, we may approximate
$\sum_{\nu}a_{\nu}a^{\dagger}_{\nu}=N$ so on the given Hilbert space,
$\mathbb{E}[H(\eta)^{2}]$ is a scalar equal to $(p/2)!{{n_{bos}}\choose
p/2}N^{p/2}+O(N^{p/2-1})$. In general, for the complex or real ensemble,
symmetrized or not, we find that $\mathbb{E}[H(\eta)^{2}]$ is a scalar equal
to $J(p/2)!{{n_{bos}}\choose p/2}N^{p/2}$, where $J$ obeys the claims of the
lemma. To verify that $\mathbb{E}[H(\eta)^{2}]$ is a scalar and to compute the
scalar to all orders in $N$, commute the annihilation operators $a_{\nu}$ to
the right. The result is some linear combination of operators with all
annihilation operators to the right of creation operators which can be written
in the form
$\sum_{\mu_{1},\ldots,\mu_{k}}(\prod_{i=1}^{k}a^{\dagger}_{\mu_{k}})(\prod_{i=1}^{k}a_{\mu_{k}})$,
and each such operator is equal to $k!{{n_{bos}}\choose k}$.

Hence, from Eq. (30),
$\partial_{\tau}^{2}\log(Z(\tau,p,N,{n_{bos}})\leq(J/2){n_{bos}}^{p/2}N^{p/2}$.
∎

Remark: this result is clearly not tight in the regime where random matrix
theory is accurate (${n_{bos}}=p/2$). It is interesting to see what happens
there. The correlation function $\mathbb{E}_{G,\eta}[{\rm
tr}\Bigl{(}\exp\\{(1-s_{1})\tau H(G)\\}H(\eta)\exp\\{s_{1}\tau
H(G)\\}H(\eta)\Bigr{)}]$ is not independent of $s_{1}$, but rather decays as a
function of $s_{1}$ for $s_{1}\leq 1/2$ (of course, it increases again as
$s_{1}$ becomes larger than $1/2$). Considering the regime $s_{1}\ll 1/2$,
using the square-root singularity at the edge of the Wigner semi-circle we can
estimate that it decays as $s_{1}^{-3/2}$. This means that the integral of
this correlation function over $s_{1}$ is dominated by its value for small
$s_{1}$ of order $1/\tau$ so that for $\tau$ large compared to the inverse
width of the semi-circle (though of course $\tau$ not too large) the integral
becomes of order $1/\tau$. This is stronger than the upper bounds here where
we have bounded by the integral by something independent of $\tau$. We may
guess that a tighter analysis will show that a similar effect will happen in
the case of ${n_{bos}}\gg p/2$; however, an important difference occurs. If we
take an eigenstate of $H(G)$ with some eigenvalue $\lambda_{0}$, and apply
$H(\eta)$ for random $\eta$, this only changes $p/2$ out of the ${n_{bos}}$
qudits in the state. So, one might guess that the resulting state will have
expectation value of $H(G)$ that is $\lambda_{0}(1-p/(2{n_{bos}}))$ rather
than an (as in the random matrix case) expectation value of $H(G)$ which is
zero. So, we may guess that the correlation function will be non-negligible
for $s_{1}\lesssim({n_{bos}}/p)\tau^{-1}$. A heuristic estimate in this
fashion suggests that the lemma below for the eigenvalue is tight up to
logarithmic factors.

From lemma 4, the following lemma is an immediate corollary:

######  Lemma 5.

Let $\lambda_{1}$ be the largest eigenvalue of $G$. Let

$E_{max}=\sqrt{2J\log(N)}{n_{bos}}^{p/4+1/2}N^{p/4}.$

Then, for any $x$,

${\rm Pr}[\lambda_{1}\geq x]\leq\exp\Bigl{(}-\frac{x-E_{max}}{\xi}\Bigr{)},$
(31)

with

$\xi=\frac{\sqrt{J}{n_{bos}}^{p/4-1/2}N^{p/4}}{\sqrt{2\log(N)}}$ (32)

So, for any $E^{\prime}$ which is $\omega(\xi)$, with high probability
$\lambda_{1}\leq E_{max}+E^{\prime}$.

###### Proof.

We have ${\rm tr}(\exp\\{\tau H(G)\\})\geq\exp(\tau\lambda_{1})$. Hence, for
any $x$, ${\rm Pr}[\lambda_{1}\geq x]\leq Z(\tau,p,N,{n_{bos}})/\exp(\tau x)$.
Since $D(N,{n_{bos}})\leq N^{n_{bos}}$, minimizing over $\tau$, we find that

${\rm Pr}[\lambda_{1}\geq
x]\leq\exp\Bigl{(}{n_{bos}}\log(N)-\frac{x^{2}}{2J{n_{bos}}^{p/2}N^{p/2}}\Bigr{)}.$

For $x=E_{max}$, the right-hand is equal to $1$ and for $x>E_{max}$ the right-
hand side decays exponentially. ∎

###  4.2 Odd $p$

We now consider the case of odd $p$. Let $Z(\tau,p,N,{n_{bos}})$ denote
$\mathbb{E}[{\rm tr}(\exp(\tau H(G)))]$ for a tensor $G$ of rank $p$, with
entries of $G$ chosen from the Gaussian distribution, with given
$N,{n_{bos}}$. In this subsection, $G$ is complex and not symmetrized. The
Hamiltonian $H(G)$ is given by Eq. (6) for even $p$ and by Eq. (3.1) for odd
$p$.

We will reduce the calculation for odd $p$ to the case for even $p$, up to a
bounded error, showing the following

######  Lemma 6.

For odd $p$, for ${n_{bos}}^{p-1}\tau N^{1/3}=o(1)$,

$Z(\sqrt{2N}\tau,2(p-1),N,{n_{bos}})\leq Z(\tau,p,N,{n_{bos}})\leq
Z(\sqrt{2N}\tau,2(p-1),N,{n_{bos}})\exp(o(N^{(p-1)/2})).$ (33)

All occurrences of $Z(\cdot)$ in the above equation refer to the complex
ensemble without symmetrizing.

Remark: the assumptions require that $\tau$ is $o(1)$; however, that is still
large enough $\tau$ to be useful later in bounding the spectrum since the
largest eigenvalues of $H(G)$ are typically large compared to $1$.

###### Proof.

For arbitrary $H(G)$, the exponential $\exp\\{\tau H(G)\\}$ can be expanded as
a series

$1+\tau H(G)+\frac{\tau^{2}}{2}H(G)^{2}+\ldots.$

We bring the expectation value (over $G$) inside the trace, and compute the
expectation value of this series. This expectation value can be computed using
Wick’s theorem to compute the expectation value of a moment of a Gaussian
distribution.

For odd $p$, each term in $H$ is a sum of two terms, one depending on the
product of two tensors $G$ and one depending on the product of two tensors
$\overline{G}$, where the overline denotes complex conjugation. Hence, the
$m$-th order term $\frac{\tau^{m}}{m!}H(G)^{m}$ is a sum of $2^{m}$ terms,
corresponding to $m$ different choices of $G$ or $\overline{G}$ in each
$H(G)$); we call each such choice a history. Each term is a product of
creation and annihilation operators depending on a product of $2m$ tensors,
some of which are $G$ and some of which are $\overline{G}$. This expectation
value of such a term is non-vanishing only if there are $m$ tensors $G$ and
$m$ tensors $\overline{G}$. In that case, the expectation value is given by
summing over ways of pairing each tensor $G$ with a distinct tensor
$\overline{G}$. There are $m!$ such pairings. Then, for each such pairing, one
computes an operator by taking the operator $\frac{\tau^{m}}{m!}H(G)^{m}$ and
replacing, for every pair, the two tensors
$G_{\mu_{1},\ldots,\mu_{p}}\overline{G}_{\nu_{1},\ldots,\nu_{p}}$ in that pair
with a product of $\delta$-functions

$\prod_{a=1}^{p}\delta_{\mu_{a},\nu_{a}}=\mathbb{E}[G_{\mu_{1},\ldots,\mu_{p}}\overline{G}_{\nu_{1},\ldots,\nu_{p}}].$

Summing this operator over pairings gives the expectation value.

Note also that here we have not symmetrized $G$. If we were to symmetrize $G$,
then
$\mathbb{E}[G_{\mu_{1},\ldots,\mu_{p}}\overline{G}_{\nu_{1},\ldots,\nu_{p}}]$
is given by $1/p!$ times a sum of $p!$ different products of
$\delta$-functions, i.e. $\sum_{\pi}\delta_{\mu_{a},\nu_{\pi(a)}}$, where
$\pi$ is a permutation. This would lead to additional terms that we need to
compute and would make the analysis more difficult (though in practice may
lead to better performance).

For given $m$, given history and given pairing, let us define a cluster as
follows: define a graph with $2m$ vertices, each vertex corresponding to a
single tensor, either $G$ or $\overline{G}$. Let there be an edge between any
two tensors which both appear in the same term in the Hamiltonian. Let there
also be an edge between any two tensors which are in a pair. Hence, this is a
graph of degree $2$ (we allow the possibility of multiple edges connecting two
vertices if we pair two tensor which both appear in the same term in the
Hamiltonian). A cluster is a connected component of this graph.

We refer to a cluster containing four vertices as a minimal cluster. A cluster
with six or more vertices is called a non-minimal cluster. Note that there are
no clusters containing only two terms because each term in $H(G)$ depends on a
product of two tensors $G$ or two tensors $\overline{G}$; if we had not taken
the complex ensemble and instead taken $G$ to be real, we would instead have
these clusters with two terms. We discuss the case of the real ensemble
further after the proof of the lemma.

The minimal clusters will turn out to give the dominant order contribution in
an expansion in $N$. The non-minimal clusters will be subleading order.

We have expressed $\frac{\tau^{m}}{m!}H(G)^{m}$ as a sum over histories and
pairings and so $\mathbb{E}[\exp(\tau H(G))]$ is a sum over $m$, histories,
and pairings. Each term in the sum over $m$, histories, and pairings for
$\mathbb{E}[\exp(\tau H(G))]$ is an operator. Hence, $Z(\tau,p,N,{n_{bos}})$
is also given by a sum over $m$, histories, and pairings as one may take the
trace of each term in $\mathbb{E}[\exp(\tau H(G))]$. Note that each $m$,
history, and pairing gives a non-negative contribution to
$Z(\tau,p,N,{n_{bos}})$, i.e., it has a non-negative trace.

Lower Bound— Now we prove the first inequality in Eq. (33), namely that
$Z(\sqrt{2N}\tau,2(p-1),N,{n_{bos}})\leq Z(\tau,p,N,{n_{bos}})$. To do this,
we first define a similar summation over $m$, histories, and pairings to
compute $Z(\tau,q,N,{n_{bos}})$ for even $q=2(p-1)$ and then we compare the
two summations.

For even $q$, each term in $H$ is a sum of two terms, one depending linearly
on a tensor $G$ and one depending on tensor $\overline{G}$. As before, the
$m$-th order term $\frac{\tau^{m}}{m!}H(G)^{m}$ is a sum of $2^{m}$ histories,
with each history corresponding to $m$ different choices of $G$ or
$\overline{G}$ in each $H(G)$. Each term is a product of creation and
annihilation operators depending on a product of $m$ tensors, some of which
are $G$ and some of which are $\overline{G}$; note the difference here from
the odd case as now there are only $m$ tensors, rather than $2m$. This
expectation value of such a term is non-vanishing only if there are $m/2$
tensors $G$ and $m/2$ tensors $\overline{G}$. In that case, the expectation
value is given as before by summing over pairings and replacing the two
tensors in the pair with $\delta$-functions.

We claim that if we consider $Z(\tau,p,N,{n_{bos}})$ for odd $p$ and consider
only the sum of pairings for which all clusters are minimal, this gives
precisely the sum of terms for $Z(\sqrt{2N}\tau,2(p-1),N,{n_{bos}})$. Since
all terms contribute non-negatively to the trace, this will prove the first
inequality.

To show this claim, let $L_{2(p-1)}$ label a given choice of $m$, history, and
pairing for $Z(\sqrt{2N}\tau,2(p-1),N,{n_{bos}})$ and let $L_{p}$ label a
given choice of $m$, history, and clusters for $Z(\tau,p,N,{n_{bos}})$ for
which all clusters are minimal. There are $2^{m/2}$ different pairings for the
given choice of clusters labelled by $L_{p}$.

We will construct a one-to-one correspondence between $L_{p}$ and $L_{2(p-1)}$
and show that the sum of the terms labelled by $L_{p}$ is equal to the term
labelled by $L_{2(p-1)}$, i.e., that they give the same operator, and hence
have the same trace. Consider given $L_{p}$. Then, $L_{2(p-1)}$ is as follows.
The value of $m$ is the same. The history is also the same: for a given
sequence of choices of $G$ or $\overline{G}$ for $L_{p}$, we use the same
sequence for $L_{2(p-1)}$. Note that in $L_{p}$, each of the $m$ choices of
$G$ or $\overline{G}$ denotes a choice that both tensors in $H(G)$ are equal
to $G$ or both are equal to $\overline{G}$, while in $L_{2(p-1)}$ each of the
$m$ choices is only a choice about a single tensor in $H(G)$.

The pairing is as follows. We introduce notation to label the tensors in
$H(G)^{m}$. For even $p$, we label the tensors by an integer
$i\in\\{1,2,\ldots,m\\}$ depending which of the $m$ factors of $H(G)$ it is
in. Here, we mean that the tensors appear in the product $H(G)^{m}=H(G)\cdot
H(G)\cdot\ldots\cdot H(G)$; we label each of the different factors in the
product $1,2,\ldots,m$ in sequence. For odd $p$, we label the tensors by a
pair $(i,w)$ for $i\in\\{1,2,\ldots,m\\}$ and $w\in\\{1,2\\}$. The integer $i$
labels which of the $m$ factors of $H(G)$ it is in and $w$ labels whether it
is the first or second of the two tensors in $H(G)$. For a pairing for
$Z(\tau,p,N,{n_{bos}})$ with all clusters minimal, then each cluster is of the
form that for some $i,j$ that we pair $(i,1)$ with $(j,1)$ and $(i,2)$ with
$(j,2)$ or of the form that we pair $(i,1)$ with $(j,2)$ and $(i,2)$ with
$(j,1)$. For a given choice of clusters, the corresponding pairing for
$L_{2(p-1)}$ pairs $i$ with $j$.

We sketch the claim about the terms. For a term in $L_{p}$, for each cluster
we have two tensors $T$ and two tensors $\overline{T}$. We replace these
tensors with the expectation value

$\mathbb{E}[\sum_{\mu_{1},\ldots,\mu_{p-1}}\sum_{\nu_{1},\ldots,\nu_{p-1}}\sum_{\sigma}T_{\mu_{1},\mu_{2},\ldots,\mu_{p-1},\sigma}T_{\nu_{1},\nu_{2},\ldots,\nu_{p-1},\sigma}\sum_{\alpha_{1},\ldots,\alpha_{p-1}}\sum_{\beta_{1},\ldots,\beta_{p-1}}\sum_{\overline{\sigma}}\overline{T}_{\alpha_{1},\alpha_{2},\ldots,\alpha_{p-1},\overline{\sigma}}\overline{T}_{\beta_{1},\beta_{2},\ldots,\beta_{p-1},\overline{\sigma}}],$

which is equal to some product of $\delta$-functions. This expectation value
then multiplies the operators
$\Bigl{(}\prod_{i=1}^{(p-1)/2}a^{\dagger}_{\mu_{i}}a^{\dagger}_{\nu_{i}}\Bigr{)}\Bigl{(}\prod_{i=(p-1)/2+1}^{p-1}a_{\mu_{i}}a_{\nu_{i}}\Bigr{)}$
and
$\Bigl{(}\prod_{i=1}^{(p-1)/2}a^{\dagger}_{\alpha_{i}}a^{\dagger}_{\beta_{i}}\Bigr{)}\Bigl{(}\prod_{i=(p-1)/2+1}^{p-1}a_{\alpha_{i}}a_{\beta_{i}}\Bigr{)}$
inserted at the appropriate places into the product $H(G)^{m}$. The
$\delta$-functions constrain $\sigma=\overline{\sigma}$; summing over this
gives a factor of $N$ for each cluster, while there are two pairings for each
cluster giving another factor of $2$ for each cluster. The number of clusters
is equal to $m/2$, giving an overall factor $(2N)^{m/2}$.

This proves the first inequality.

Upper bound— Now we prove the second inequality in Eq. (33). To do this, we
define the following quantity for $q$ which is a multiple of $4$ (note that
$q=2(p-1)$ is a multiple of $4$ if $p$ is odd):

$Z^{\prime}(\tau,\tau^{\prime},q,N,{n_{bos}})\equiv\mathbb{E}[{\rm
tr}(\exp\\{\tau H(G)+\tau^{\prime}H(G^{\prime})^{2}\\})],$ (34)

where the expectation value is over tensors $G$ of order $q$ chosen from the
complex ensemble and tensors $G^{\prime}$ of order $q/2$ chosen also from the
complex ensemble (note that $q/2$ is even). Note that we square
$H(G^{\prime})$ in the exponent.

We will prove that

$Z(\tau,p,N,{n_{bos}})\leq
Z^{\prime}(\sqrt{2N}\tau,\tau^{\prime},2(p-1),N,{n_{bos}})$ (35)

for $\tau^{\prime}=N^{1/3}\tau$. From this, the second inequality in Eq. (33)
follows. To see this, we have, for any $G^{\prime}$,
$\|H(G^{\prime})^{2}\|\leq O({n_{bos}}^{p-1})\|G^{\prime}\|^{2}$ where the
operator norm $\|H(G^{\prime})^{2}\|$ denotes the largest eigenvalue in
absolute value and $\|G^{\prime}\|$ denotes the largest singular value of
$G^{\prime}$, regarding $G^{\prime}$ as a matrix of size
$N^{q/4}$-by-$N^{q/4}$. So, by the Golden-Thompson inequality,

$\displaystyle\mathbb{E}_{G}[{\rm tr}(\exp\\{\tau
H(G)+\tau^{\prime}H(G^{\prime})^{2}\\})]$ $\displaystyle\leq$
$\displaystyle\mathbb{E}_{G}[{\rm tr}(\exp\\{\tau
H(G)\\})]\mathbb{E}_{G^{\prime}}[\exp(O({n_{bos}}^{p-1})\tau^{\prime}\|G^{\prime}\|^{2})]$
$\displaystyle=$ $\displaystyle
Z(\sqrt{2N}\tau,2(p-1),N,{n_{bos}})\mathbb{E}_{G^{\prime}}[\exp(O({n_{bos}}^{p-1})\tau^{\prime}\|G^{\prime}\|^{2})].$

where the subscript $G$ or $G^{\prime}$ denotes the expectation over $G$ or
$G^{\prime}$. The matrix $G^{\prime}$ has complex entries; we can bound its
operator norm by the sum of the operator norms of its Hermitian and anti-
Hermitian parts. For ${n_{bos}}^{p-1}\tau^{\prime}=o(1)$, we have
$\mathbb{E}_{G^{\prime}}[\exp(O({n_{bos}}^{p-1})\tau^{\prime}\|G^{\prime}\|^{2})]=\exp[O({n_{bos}}^{p-1})\tau^{\prime}O(N^{q/4})]$
as can be shown using Hermite polynomials [12] to compute the probability
distribution of the eigenvalues of a random matrix and using the decay of a
Hermite polynomial multiplying a Gaussian (the study of the largest eigenvalue
of a random matrix is a rich field if we consider the probability distribution
near the edge of the Wigner semi-circle but here we are satisfied to consider
the probability distribution for eigenvalues which are some constant factor
$>1$ multiplying the edge of the semi-circle).

To show Eq. (35), let $L_{p}$ label the set of terms for
$Z(\tau,p,N,{n_{bos}})$ with a given choice of $m$, history, clusters (the
clusters need not be minimal), and choice of pairing for all tensors in non-
minimal clusters (we do not specify the pairing of the tensors in the minimal
clusters; if there are $n_{min}$ minimal clusters then there are $2^{n_{min}}$
terms in the set). Let $L^{\prime}_{2(p-1)}$ label a term for
$Z^{\prime}(\sqrt{2N}\tau,\tau^{\prime},2(p-1),N,{n_{bos}})$ with given $m$,
history, clusters, and pairing. Here, for $Z^{\prime}$, a history corresponds
to $m$ choices of either $H(G)$ or $H^{\prime}(G)^{2}$ and also to choices of
$G$ or $\overline{G}$, or $G^{\prime}$ or $\overline{G}^{\prime}$, for each
$H(G)$ or $H(G^{\prime})$. We define a map from choice of $L_{p}$ to choice of
$L^{\prime}_{2(p-1)}$ so that the sum of terms labelled by $L_{p}$ is equal to
the term labelled by the corresponding $L^{\prime}_{2(p-1)}$. This map will be
one-to-one but it will not be onto. However, since all terms are non-negative,
this will establish the needed inequality.

The map is as follows. The $m$ will be the same. Using the notation above, if
tensor $(i,1)$ is in a minimal cluster (and hence, so is $(i,2)$) in the set
labelled by $L_{p}$, then in the history for $L^{\prime}_{2(p-1)}$ for the
$i$-th term we choose $\tau H(G)$, while if $(i,1)$ is not in a minimal
cluster, then we choose $\tau^{\prime}H(G^{\prime})^{2}$. Given a history
labelled by $L_{p}$, corresponding to a choice of tensor or its complex
conjugate (i.e., either $G$ or $\overline{G}$) in each of the $m$ different
$H(G)$ in a term in the series for $Z(\tau,p,N,{n_{bos}})$, then we make the
same choice of tensor or its complex conjugate in each of the $m$ different
$H(G)$ or $H(G^{\prime})^{2}$ in the history labelled by
$L^{\prime}_{2(p-1)}$. That is, if we choose $G$ in the $i$-th $H(G)$ in some
terms in the series for $Z(\tau,p,N,{n_{bos}})$, then we choose $G$ in $H(G)$
or $G^{\prime}$ in both terms in $H(G^{\prime})^{2}$ and similarly if we
choose $\overline{G}$ we choose $\overline{G}$ in $H(G)$ and
$\overline{G}^{\prime}$ in both terms in $H(G^{\prime})^{2}$.

We label the tensors in the expansion for
$Z^{\prime}(\sqrt{2N}\tau,\tau^{\prime},2(p-1),N,{n_{bos}})$ either by an
integer $i$, if it appears in $H(G)$, or by a pair $(i,w)$, if it appears in
$H(G^{\prime})^{2}$, in which case the index $w\in\\{1,2\\}$ labels which of
the two $H(G^{\prime})$ it is in. Finally, we define the pairing labelled by
$L^{\prime}_{2(p-1)}$. For a minimal cluster labelled by $L_{p}$, pairing
$(j,1)$ and $(i,2)$ with $(j,2)$ or $(i,1)$ with $(j,2)$ and $(i,2)$ with
$(j,1)$, then in the pairing labelled by $L^{\prime}_{2(p-1)}$ we pair $i$
with $j$..

If a cluster is non-minimal, then we simply use the same pairing for the
corresponding tensors in $L^{\prime}_{2(p-1)}$. That is, suppose a cluster
pairs $(i_{1},w_{1})$ with $(i_{2},w^{\prime}_{2})$, and pairs $(i_{2},w_{2})$
with $(i_{3},w^{\prime}_{3})$, and so on, where $w^{\prime}_{a}=1$ if
$w_{a}=2$ and $w^{\prime}_{a}=2$ if $w_{a}=1$. Then, we also pair
$(i_{1},w_{1})$ with $(i_{2},w^{\prime}_{2})$, and pairs $(i_{2},w_{2})$ with
$(i_{3},w^{\prime}_{3})$, and so on.

The smallest non-minimal cluster has six vertices. In every cluster, minimal
or not, there is a sum over some index (for example, the sum over the index
$\sigma=\sigma^{\prime}$ in the lower bound calculation above) which gives a
factor $N$. Thus, taking $\tau^{\prime}=\tau N^{1/3}$ accounts for this
factor. No factor of $2$ occurs for the non-minimal clusters. ∎

Remark: since we have chosen the entries of $G$ from the complex ensemble, we
have that $\mathbb{E}[H(G)]=0$ for $p$ odd. If instead we had chosen the
entries of $G$ from the real ensemble we would have (considering the specific
case $p=3$ for simplicity and not symmetrizing $G$, again for simplicity) a
non-vanishing expectation value since

$\mathbb{E}[T_{\mu_{1},\mu_{2},\sigma}T_{\nu_{1},\nu_{2},\sigma}]=N\delta_{\mu_{1},\nu_{1}}\delta_{\mu_{2},\nu_{2}},$
(37)

so that $\mathbb{E}[H(G)]\propto
N\sum_{\mu_{1},\mu_{2}}a^{\dagger}_{\mu_{2}}a^{\dagger}_{\mu_{2}}a_{\mu_{1}}a_{\mu_{1}}$.
Such a term is sometimes called a pairing term or a “cooperon” in the study of
disordered system in physics. In the case ${n_{bos}}=2$ (the smallest possible
for $p=3$), this term has operator norm $N^{2}$. This is much larger than
$N^{1/2}$ times the expected operator norm of $H(G)$ for $p=2(p-1)=4$, i.e.,
that expected operator norm is proportional to $N$ by random matrix theory,
and $N^{2}\gg N^{3/2}$.

There may be other ways to deal with this non-vanishing expectation value
other than using the complex ensemble. One way is to use the real ensemble,
but to consider the Hamiltonian $H(T_{0})-M$, where we define
$M=N\sum_{\mu_{1},\mu_{2}}a^{\dagger}_{\mu_{2}}a^{\dagger}_{\mu_{2}}a_{\mu_{1}}a_{\mu_{1}}$
for $p=3$. In this case, the added term $-M$ cancels the expectation value of
$H(G)$ term-by-term in the perturbation expansion. However, if we do this we
still have some additional terms when we consider clusters of four or more
vertices. We expect that the clusters of six or more vertices are still
negligible, but the structure of the clusters of four vertices becomes more
complicated. We leave the analysis of this case for the future, but we expect
that it would work and may be practically useful.

From lemma 6, the following lemma is an immediate corollary:

######  Lemma 7.

Let $\lambda_{1}$ be the largest eigenvalue of $G$. Let

$E_{max}=2\sqrt{J\log(N)}{n_{bos}}^{p/2}N^{p/2},$

where $J$ is the $J$ of lemma 4 for $Z(\tau,2(p-1),N,{n_{bos}})$. Then, for
any $x$,

${\rm Pr}[\lambda_{1}\geq x]\leq\exp\Bigl{(}-\frac{x-E_{max}}{\xi}\Bigr{)},$
(38)

with

$\xi=\frac{\sqrt{J}{n_{bos}}^{p/2-1}N^{p/2}}{\sqrt{\log(N)}}.$ (39)

So, for any $E^{\prime}$ which is $\omega(\xi)$, with high probability
$\lambda_{1}\leq E_{max}+E^{\prime}$.

###### Proof.

From lemmas 4,6, for ${n_{bos}}^{p-1}\tau N^{1/3}=o(1)$, we have

$Z(\tau,p,N,{n_{bos}})\leq\exp(\tau^{2}J{n_{bos}}^{p-1}N^{p})\exp(o(N^{(p-1)/2})).$

Let $\tau={n_{bos}}^{1-p/2}N^{-p/2}\sqrt{\log(N)}/\sqrt{J}$. For
${n_{bos}}=o(N^{(p/2-1/3)/(1+p/2)}\log(N)^{1/(2+p)}$, the condition
${n_{bos}}^{p-1}\tau N^{1/3}=o(1)$ holds. So, after some algebra, for any $x$,
the result follows. ∎

It is worth noting that lemma 6 has the following corollary:

######  Corollary 2.

For odd $p$, for ${n_{bos}}=p-1$, with high probability
$\lambda_{1}=O(N^{p/2})$.

###### Proof.

To prove this we use the existence of tighter bounds for even $q=2(p-1)$ when
${n_{bos}}=q/2$. By Eq. (33), for ${n_{bos}}^{p}\tau N^{1/3}=o(1)$, we have
$Z(\tau,p,N,{n_{bos}})\leq
Z(\sqrt{2N}\tau,2(p-1),N,{n_{bos}})\exp(o(N^{(p-1)/2}))$. Since we are
considering fixed ${n_{bos}}$, this holds for $\tau N^{1/3}=o(1)$. We have
that $Z(\sqrt{2N}\tau,2(p-1),N,p-1)=\mathbb{E}_{G}[\exp(\sqrt{2N}\tau H(G)]$,
but the Hamiltonian $H(G)$ is a random matrix of size $N^{p-1}$-by-$N^{p-1}$
chosen from the so-called Gaussian unitary ensemble [12]. With high
probability, the largest eigenvalue of this matrix is $\Theta(N^{(p-1)/2})$.
Further, we can choose $\tau=\omega(\log(N)N^{-p/2})$ so that
$Z(\sqrt{2N}\tau,2(p-1),N,p-1)=N\exp(O(\tau N^{p/2}))$; for example, this can
be shown by using orthogonal polynomials to bound the probability distribution
of the largest eigenvalue as discussed previously. We have ${\rm
Pr}[\lambda_{1}>x]\leq Z(\tau,p,N,{n_{bos}})\exp(-\tau x)$. For $x$
sufficiently large compared to $N^{p/2}$, for $\tau=\omega(\log(N)N^{-p/2})$,
the right-hand side is $o(1)$. ∎

##  5 Quantum and Classical Algorithms

We now discuss the complexity of classical and quantum algorithms to implement
the needed linear algebra. In particular, we need to determine if that largest
eigenvalue is larger than $E_{cut}$ and we need to find some vector in the
eigenspace of eigenvalue $\geq E_{cut}$. We emphasize that it is not necessary
to find the leading eigenvector itself.

We will use $\psi_{\rm target}$ to denote this leading eigenvector. Note that
if we lower bound the squared overlap of some vector with $\psi_{\rm target}$,
this will lower bound the probability of success of phase estimation.

We have defined algorithm 2 so that in the detection step it uses a hard
cutoff on eigenvalue: if the leading eigenvalue is $\geq E_{cut}$ it reports
that $\lambda=\overline{\lambda}$ while otherwise it reports that $\lambda=0$.
However, no algorithm, classical or quantum, will be able to compute the
leading eigenvalue exactly; there will always be some limit on the precision.
Fortunately, the proofs of theorems 2,4 show that if $E_{0}\geq(1+c)E_{max}$
for any $c>0$, then if $\lambda=\overline{\lambda}$ then with high probability
$\lambda_{1}\geq(1-\eta)E_{0}+\eta E_{max}$ for any $0<\eta<1$ while if
$\lambda=0$ then with high probability
$\lambda_{1}\leq(1-\eta^{\prime})E_{0}+\eta^{\prime}E_{max}$ for any
$0<\eta^{\prime}<1$. For example, we might take $\eta=1/8$ and
$\eta^{\prime}=1/2$. So, it suffices instead to implement some “soft” estimate
of the leading eigenvalue which will be very likely to give one result (i.e.,
reporting that $\lambda=\overline{\lambda}$) if the leading eigenvalue is
larger than $(7/8)E_{0}+(1/8)E_{max}$ but very unlikely to give that result if
the leading eigenvalue is $\leq E_{cut}$.

In the quantum algorithms, to obtain some vector in the eigenspace of
eigenvalue $\geq E_{cut}$ and to do this soft estimate, we will implement an
approximate projector by phase estimation in the quantum algorithms onto the
eigenspace of eigenvalue $\geq(E_{0}+E_{cut})/2=(3/4)E_{0}+(1/4)E_{cut}$. By
doing this, the phase estimation error will become negligible when considering
the projection of the resulting vector onto the eigenspace with eigenvalue
$<E_{cut}$. Similarly, in the classical power method we will take the number
of iterations sufficiently large that the vector has negligible projection on
the eigenspace with eigenvalue $<E_{cut}$ and further so that it has
expectation value for $H(T_{0})$ greater than $E_{cut}$.

We begin with a description of some classical algorithms. The time and space
requirements for the classical algorithms are of course not intended to
represent a lower bound; rather, they represent times that can be achieved
using standard algorithms in the literature. We then give quantum algorithms.
Finally, we give a further improvement to the quantum algorithm that may be
useful in practice.

When we refer to “space” in a classical algorithm, if we store a
$D(N,{n_{bos}})$-dimensional vector, the space requirement is equal to
$D(N,{n_{bos}})$ multiplied by the number of bits to store a single entry of
the vector. In the classical algorithms, we will not discuss issues with
finite precision arithmetic in detail. Since we will be applying operators of
the form $H(T_{0})^{m}$ to vectors in the “path integral” methods, we might
need to approximate each entry of the vector to accuracy
$D(N,{n_{bos}})^{-1}\|H(T_{0})\|^{-m}=O({\rm poly}(N^{{n_{bos}}m}))^{-1}$.
However, the required number of bits is then only $O(m{n_{bos}}\log(N))$ and
$m$ will be logarithmic in $N$ so the required number of bits will be only
polylogarithmic in $N$.

###  5.1 Classical Algorithms

Classically, the most obvious algorithm is to perform an eigendecomposition on
$H(T_{0})$. This requires storing matrices of size
$D(N,{n_{bos}})$-by-$D(N,{n_{bos}})$ so that the space required is
$\tilde{O}(D(N,{n_{bos}})^{2})$ and the time is $D(N,{n_{bos}})^{\omega}$
where $\omega$ is the matrix multiplication exponent [13], though of course in
practice the time is closer to $D(N,{n_{bos}})^{3}$.

However, there is no need to perform a full eigendecomposition. One can
instead initialize a random vector and then apply the power method to extract
some eigenvector of $H(T_{0})$ in the eigenspace with eigenvalue $\geq
E_{cut}$. The space required is then only $\tilde{O}(D(N,{n_{bos}}))$. The
time required for a single iteration of the power method is
$\tilde{O}(D(N,{n_{bos}}))$. If $\lambda=\overline{\lambda}$, then in
$O(\log(D(N,{n_{bos}}))/\log(E_{0}/E_{max}))$ iterations, the resulting vector
will have an $1-o(1)$ projection onto the eigenspace with eigenvalue
$\geq(E_{0}+E_{cut})/2$ and a negligible projection onto the eigenspace with
eigenvalue $<E_{cut}$. So, after this many iterations, one can compute the
expectation value of $H(T_{0})$ on that vector to perform detection, i.e., if
$\lambda=\overline{\lambda}$ the expectation will be larger than $E_{cut}$ but
if $\lambda=0$ the expectation will be close to $E_{max}$, and one can compute
the single particle density matrix of the vector after those iterations. So,
the time is $\tilde{O}(D(N,{n_{bos}}))O(1/\log(E_{0}/E_{max}))$.

This power method still requires exponential (in ${n_{bos}}$) space. In fact,
one can use only polynomial space. One obvious choice is to perform a “path
integral”. That is, given a random initial state $\Psi_{\rm random}$, the
power method computes the state $H(T_{0})^{m}|\Psi_{\rm random}\rangle$, for
some exponent $m$ which gives the number of iterations in the power method.
So, we wish to compute

$\frac{\langle\Psi_{\rm
random}|H(T_{0})^{m}a^{\dagger}_{\mu}a_{\nu}H(T_{0})^{m}|\Psi_{\rm
random}\rangle}{\langle\Psi_{\rm random}|H(T_{0})^{2m}|\Psi_{\rm
random}\rangle},$

where the denominator is to correctly normalize the state. Let us choose the
initial state to be a (normalized) random state from the basis for the
symmetric basis given before. We make this choice to make the “path integral”
simpler and since this is a complete orthonormal basis of states, a random
state from this basis has expected overlap with the largest eigenvector of
$H(T_{0})$ equal to $1/N^{{n_{bos}}}$. Then, both the numerator and
denominator above can be expressed by a summation over intermediate states
from this basis, requiring only space $\tilde{O}(\log(D(N,{n_{bos}}))m)$; this
summation is a “path integral”. The time required however is now
$\tilde{O}(D(N,{n_{bos}})^{m})$ and so becomes significantly worse if the
number of iterations is much more than $1$.

An improvement to this time can be obtained by using the algorithm of theorem
4.1 of Ref. [11]. This algorithm is expressed in terms of qubits, but we can
translate the algorithm to get an algorithm for a single qudit of dimension
$D(N,{n_{bos}})$, or, using the tradeoffs discussed there, ${n_{bos}}$ qudits
each of dimension $N$.

Let us explain this in detail, following the discussion there. We wish to
compute $\langle y|C|x\rangle$, for some fixed basis states $x,y$ from the
basis above (in our application $x=y$, where $C$ is a “circuit” of some depth
$d$ (for example $C=2m+1$ for the numerator above and $C=m$ for the
denominator above). While in Ref. [11], a circuit was assumed to be unitary,
there in fact is no need to assume that; for non-unitary circuits, the issues
with finite precision arithmetic do become more important but as remarked
above this can be handled with only polylogarithmic overhead. For us, a
circuit is not necessarily unitary; rather it is built out of a sequence of
operations, each of which is either $H(T_{0})$ or $a^{\dagger}_{\mu}a_{\nu}$;
more generally a circuit might include any operations built out of some tensor
of order $O(1)$ multiplying some number of creation and annihilation
operators.

For $d=1$, $\langle y|C|x\rangle$ can be computed computed in time ${\rm
poly}({n_{bos}})$. For $d>1$, we have

$\displaystyle\langle y|C|x\rangle=\sum_{z}\langle y|C_{[d\leftarrow
d/2+1]}|z\rangle\cdot\langle z|C_{[d/2\leftarrow 1]}|x\rangle,$ (40)

where the summation is over states $z$ in the basis and where $C_{[d\leftarrow
d/2+1]}$ and $C_{[d/2\leftarrow 1]}$ are subcircuits from the second and first
half of $C$. If $F(d)$ is the runtime at depth $d$, we have $F(d)\leq 2\cdot
D(N,{n_{bos}})\cdot F(\lceil d/2\rceil)$. So, $F(d)\leq{\rm
poly}({n_{bos}})(2D(N,{n_{bos}}))^{\lceil\log(d)\rceil}.$ So, the runtime is
$\tilde{O}((2D(N,{n_{bos}}))^{\lceil\log(2m+1)\rceil})$. This is much faster
than the path integral method but potentially much slower than the power
method, depending on the required $m$, i.e., for $E_{0}/E_{max}=\Theta(1)$, we
need $m$ proportional to $\log(D(N,{n_{bos}}))$ and so the time required is
superpolynomially worse than the time required if one stores the full vector.

###  5.2 Quantum Algorithms

We now discuss quantum algorithms for the same problem. In contrast to the
classical algorithms above, all these algorithms take only polynomial space.
First let us describe a simple algorithm, given as algorithm 3. We then
describe a sequence of improvements.

For algorithm 3 and all subsequent algorithms, we analyze assuming that
$\lambda=\overline{\lambda}$, i.e., for purposes of analysis we consider the
problem of recovery rather than detection. All these algorithms report success
or failure and, if they fail, they are rerun until they succeed. We give
bounds on the expected runtime under the assumption that
$\lambda=\overline{\lambda}$. If we consider the problem of detection, and if
$\lambda=0$, then the algorithm will not report success in the given runtime
(since it will be unable to succeed in a certain phase estimation procedure)
and so all these algorithms can also be used for detection by running them for
some multiple of the given runtime and reporting that $\lambda=0$ if success
is not reported in that time.

####  5.2.1 Maximally Entangled (or Maximally Mixed) Input State

Algorithms 3,4 in this subsubsection work on a Hilbert space which is the
tensor product of two Hilbert spaces, each of which have dimension
$D(N,{n_{bos}})$. (In some versions of the Hamiltonian simulation used in the
algorithms, it is convenient to embed the symmetric subspace of dimension
$D(N,{n_{bos}})$ within the full Hilbert space.)

We use a tensor product notation $A\otimes B$ to denote an operator that is a
tensor product of two operators $A,B$ on the two different tensor factors.

Algorithm 3 Quantum Algorithm (simplest, unamplified version). This and all
other quantum algorithms have the same inputs. outputs, and parameter choices
as algorithm 2.

  * 1.

Prepare a maximally entangled state between the two qudits.

  * 2.

Apply phase estimation using $H(T_{0})\otimes I$. Let $\psi$ be the resulting
state. If the resulting eigenvalue is larger than $(E_{0}+E_{cut})/2$, report
“success”. Otherwise, report “failure”.

  * 3.

If success is reported, measure and return
$\langle\psi|a^{\dagger}_{\mu}a_{\nu}\otimes I|\psi\rangle.$

Steps $1-2$ are designed to prepare a state whose density matrix on the first
qudit has large projection onto the eigenspace of eigenvalue $\geq E_{cut}$.
For purposes of analysis, we trace out the second qudit, so that the input
state on the first qudit is a maximally mixed state. If success is reported
then (ignoring phase estimation error) we have indeed projected onto this
eigenspace. Ignoring phase estimation error, the probability of success is
$\geq 1/D(N,{n_{bos}})$.

Remark: we prepare a maximally entangled state between the two qudits so that
the density matrix on the first state is maximally mixed; we could equally
well modify the algorithm to use only a single qudit (reducing the space by a
factor of $2$) and prepare a random state on the first qudit. This
modification however requires some care when we discuss a version of the
algorithm that uses amplitude amplification below. In the unamplified version,
we can choose a new random state (for example, choosing it uniformly from any
orthogonal basis) each time we perform phase estimation; however, in the
amplified version, we must choose a fixed random initial state on the first
qudit and amplify the algorithm with that choice of initial state. We claim
(we omit the proof) that if one picks a tensor product of random single qudit
states $|v_{1}\otimes v_{2}\otimes\ldots\otimes v_{{n_{bos}}}\rangle$ where
$v_{a}$ are independently chosen from a Haar uniform distribution on the
sphere $|v_{a}|=1$, with high probability the state has squared overlap with
the leading eigenvector close to $N^{-{n_{bos}}}$; more precisely, “close”
means that with high probability the logarithm of the squared overlap is at
least $-{n_{bos}}\log(N)\cdot(1+o(1))$. Note that we might not have this
property of the overlap if we had chosen the initial state from the
computational basis, for example. Note also that we can efficiently prepare
states from this distribution. Sketch of proof of claim: consider the
logarithm of the probability that the sequence of measurements
$|v_{i}\rangle_{i}\langle v_{i}|$ succeeds for $i=1,\ldots,{n_{bos}}$ in
sequence. The probability that the $i$-th measurement succeeds, conditioned on
the previous measurements succeeding, can be computed from the trace of
$|v_{i}\rangle\langle v_{i}|$ with the reduced density matrix on the $i$-th
qudit of some state (i.e., the leading eigenvector projected by the previous
measurements), and for any such reduced density matrix with high probability
the logarithm of the trace is at least $-\log(N)\cdot(1+o(1))$.

Step $3$ of the algorithm measures some property of a state in the eigenspace
with eigenvalue $\geq E_{cut}$. It is possible that each time the algorithm is
run, one obtains a different energy measurement and hence a different state,
so that measuring this property of the state gives some expectation value of
$a^{\dagger}_{\mu}a_{\nu}$ in a mixed state. This does not matter since
theorems 3 or 5 also hold for mixed states.

We explain the measurement in step 3 in more detail below. The simplest
possibility is to simply measure one matrix element of
$a^{\dagger}_{\mu}a_{\nu}$. Since there are $N^{2}$ matrix elements, we then
need to repeat the algorithm ${\rm poly}(N,\log(\epsilon))$ times to measure
each matrix element to accuracy $\epsilon$. We explain improvements to this in
subsection 5.3.

We explain the phase estimation in more detail below. First, let us analyze
the algorithm in a rough outline. Let the phase estimation be carried out to a
precision sufficiently smaller than $E_{0}-E_{max}$. To define this, we work
in an eigenbasis of $H(T_{0})$. Let $\tilde{\epsilon}$ be a bound on the error
probability of the phase estimation. More precisely, we will say that phase
estimation implements an operator $E_{PE}$ which is diagonal in the eigenbasis
such that on a normalized eigenvector $v_{i}$ with eigenvalue $\lambda_{i}$ we
have

$\displaystyle\lambda_{i}<E_{cut}\;\rightarrow\;\langle
v_{i}|E_{PE}|v_{i}\rangle\leq\tilde{\epsilon},$ (41)
$\displaystyle\lambda_{i}>(7/8)E_{0}+(1/8)E_{max}\;\rightarrow\;\langle
v_{i}|E_{PE}|v_{i}\rangle\geq 1-\tilde{\epsilon}.$

The reader should appreciate that there are a few energy scales here, chosen
rather arbitrarily. The exact value of the energies are not too important. We
have picked $E_{cut}$ to be partway between $E_{0}$ and $E_{max}$ so that it
is very unlikely that the largest eigenvalue of $H(G)$ is above $E_{cut}$ and
also very unlikely that the $\lambda_{1}<E_{cut}$. We have picked the energy
cutoff in step 2 to be $(E_{0}+E_{cut})/2$ simply to pick some energy partway
between $E_{0}$ and $E_{cut}$ so that it is very unlikely that phase
estimation reports success on an eigenstate with energy $<E_{cut}$; see first
line of Eq. (41). In the last line of Eq. (41) we wrote
$(7/8)E_{0}+(1/8)E_{max}$ simply to pick some energy scale slightly below
$E_{0}$ above which it is very likely that phase estimation reports success;
for us later, it would suffices to have any instead a bound for
$\lambda_{i}\geq(1-\eta)E_{0}+\eta E_{max}$ for any $\eta>0$. Of course, the
two lines in Eq. (41) are not completely symmetric about $(E_{0}+E_{cut})/2$.

Then, choose $\tilde{\epsilon}=\epsilon/D(N,{n_{bos}})$, so that the algorithm
reports success with probability at least
$(1-\tilde{\epsilon})/D(N,{n_{bos}})$ and, given that the algorithm reports
success, the resulting state has projection onto the eigenspace with
eigenvalue $\geq E_{cut}$ which is greater than or equal to

$\frac{(1-\tilde{\epsilon})}{(1-\tilde{\epsilon})+(D(N,{n_{bos}})-1)\tilde{\epsilon}}=1-O(\epsilon),$

where the first term in the denominator is the probability that it reports
success on $\psi_{\rm target}$ as input and the second term is the probability
of reporting success on a state in the eigenspace with eigenvalue $<E_{cut}$,
multiplied by $D(N,{n_{bos}})-1$, i.e., multiplied by an upper bound on the
dimensionality of that eigenspace. Taking $\epsilon\ll 1$, we can obtain a
large projection onto the eigenspace with eigenvalue $\geq E_{cut}$, so that
the cost of phase estimation increases logarithmically with
$D(N,{n_{bos}})/\epsilon$.

The success probability for $\epsilon\ll 1$ is greater than or equal to
$(1/D(N,{n_{bos}}))(1-\epsilon/D(N,{n_{bos}}))$, so for small $\epsilon$ it is
very close to $1/D(N,{n_{bos}})$. Hence, repeating the algorithm until it
succeeds, the expected runtime to obtain a single measurement of one matrix
element of $a^{\dagger}_{\mu}a_{\nu}$ is bounded the time for phase estimation
multiplied by $O(D(N,{n_{bos}}))$.

To perform phase estimation, we use controlled simulation of the Hamiltonian
$H(T_{0})$. There are a large number of quantum simulation algorithms which
would work here, such as Refs. [14, 15, 16, 17, 18] to name just a few. There
are two broad possibilities. The first possibility is to work in the symmetric
subspace of dimension $D(N,{n_{bos}})$. In this case, $H(T_{0})$ is a sparse
Hamiltonian, and sparse simulation algorithms apply. The second possibility is
to use the Hamiltonian of Eq. (5) and embed the symmetric subspace into the
full Hilbert space; in this case, $H(T_{0})$ is a local Hamiltonian, in that
each terms acts on a small number of qudits, each of dimension $N$, and local
simulation algorithms apply. The cost for these algorithms to simulate for a
time $t$ to error $\tilde{\epsilon}$ is ${\rm
poly}(t\|H(T_{0})\|,{n_{bos}},N,\log(\tilde{\epsilon}))$.

Using the simplest phase estimation algorithm of Ref. [19], the number of bits
that we need to phase estimate is $s=O(\log(\|H(T_{0})\|/(E_{0}-E_{max}))$.
The most expensive bit to obtain is the least significant bit, since obtaining
the $j$-th least significant bit requires simulating for a time proportional
to $2^{s-j}(E_{0}-E_{max})^{-1}$. So, we can obtain the least significant bit
to error $\tilde{\epsilon}/2$, then obtain the next least significant bit to
error $\tilde{\epsilon}/4$, and so on, making the total error
$\tilde{\epsilon}$. Of course, a large number of variations of the Kitaev
phase estimation algorithm exist in the literature, and any could be used
here.

With high probability, $\|H(T_{0})\|$ is ${\rm poly}(N)$. Thus, with high
probability the time for phase estimation is ${\rm
poly}(N,{n_{bos}},1/(E_{0}-E_{max}),\log(D(N,{n_{bos}})/\epsilon))$, giving an
algorithm runtime

$D(N,{n_{bos}}){\rm
poly}(N,{n_{bos}},1/(E_{0}-E_{max}),\log(D(N,{n_{bos}})/\epsilon)).$

We can speed this algorithm up quadratically by applying amplitude
amplification [20]. Modify the phase estimation step $2$ of algorithm 3 so
that the algorithm phase estimates the eigenvalue, determines if the
eigenvalue is larger than $E_{cut}$, then uncomputes the eigenvalue, returning
just a single bit of success or failure. See algorithm 4. Then, applying
amplitude amplification, with high probability the algorithm succeeds in
expected time $D(N,{n_{bos}})^{1/2}{\rm
poly}(N,{n_{bos}},1/(E_{0}-E_{max}),\log(D(N,{n_{bos}})/\epsilon)).$
Multiplying by ${\rm poly}(N,1/\epsilon)$ to measure
$a^{\dagger}_{\mu}a_{\nu}$ to accuracy $\epsilon$, the the expected time is
still

$D(N,{n_{bos}})^{1/2}{\rm
poly}(N,{n_{bos}},1/(E_{0}-E_{max}),\log(D(N,{n_{bos}})/\epsilon)),$

giving a quadratic time improvement, up to ${\rm poly}(N)$ factors, and an
exponential space improvement, over the fastest classical algorithm described
above.

Algorithm 4 Quantum Algorithm (amplified version)

  * 1.

Apply amplitude amplification to steps $1-2$ of algorithm 3, modifying step
$2$ to uncompute the eigenvalue and return only success or failure.

  * 2.

If success is reported, measure and return
$\langle\psi|a^{\dagger}_{\mu}a_{\nu}\otimes I|\psi\rangle.$

####  5.2.2 Chosen Input State: Simple Version

We can obtain a further quadratic speedup by modifying the initial state that
we phase estimate. In this subsubsection, let us first explain an algorithm
that gives the basic idea of the initial state preparation; we will only be
able to prove some slightly weaker results for this algorithm (in particular
we will prove a lower bound on the average inverse runtime, rather than an
upper bound on the average runtime). We expect that this is primarily a
technical issue and that some concentration of measure argument should allow
us to prove an upper bound on the average runtime. We will then describe in
the next subsubsection a modification to the algorithm which avoids this
technical difficulty and for which we can prove a quadratic improvement
without further assumption.

Instead of a maximally entangled state, we can use the tensor $T_{0}$ to
prepare a state with a larger projection onto $\psi_{\rm target}$. In this
subsubsection, we will work in the $D^{n_{bos}}$-dimensional Hilbert space of
Eq. (5), so we will use ${n_{bos}}$ qudits each of dimension $N$. The
unamplified version is algorithm 5 and a version with amplitude amplification
is algorithm 6.

The most important new step that must be explained is the initial state
preparation (we discuss some other details at the end of this subsubsection).
We use the fact that, given a classical list of amplitudes for some
$M$-dimensional vector, with the vector having unit norm, we can prepare a
quantum state on an $M$-dimensional qudit with the given amplitudes (up to an
ill-defined overall phase, of course) using a quantum circuit of depth $O(M)$
and using $O(M)$ classical computation. For example, labelling the basis
states $|0\rangle,|1\rangle,\ldots,|M-1\rangle$, one can start with initial
state $|0\rangle$ and apply a sequence of $M-1$ rotations in the two
dimensional subspaces spanned by $i\rangle,|i+1\rangle$ for $i=0,\ldots,M-2$.

Algorithm 5 Quantum Algorithm (improved input state, unamplified version)

  * 1.

Use $T_{0}$ to prepare the initial state $\Psi_{\rm input}$ of Eq. (42).

  * 2.

If the initial state is not in the symmetric subspace, report “failure”. If
the state is in the symmetric subspace, apply phase estimation using
$H(T_{0})$. Let $\psi$ be the resulting state. If the resulting eigenvalue is
larger than $(E_{0}+E_{cut})/2$, report “success”. Otherwise, report
“failure”.

  * 3.

If success is reported, measure and return
$\langle\psi|a^{\dagger}_{\mu}a_{\nu}|\psi\rangle.$

Algorithm 6 Quantum Algorithm (amplified version)

  * 1.

Apply amplitude amplification to steps $1-2$ of algorithm 5, modifying step
$2$ to uncompute the eigenvalue and uncompute the determination of whether the
state is in the symmetric subspace and to return only success or failure.

  * 2.

If success is reported, measure and return
$\langle\psi|a^{\dagger}_{\mu}a_{\nu}|\psi\rangle.$

In this subsubsection, for simplicity in analysis, we assume that the error
$\epsilon$ in phase estimation is small enough to be negligible.

We use the same method to produce the input state for both even and odd $p$.
First consider the case that ${n_{bos}}$ is an integer multiple of $p$. For
any tensor $T$ of order $p$, let $|T\rangle$ denote the vector on $p$ qudits
(each of dimension $N$) with amplitudes given by the entries of the tensor.
This vector is correctly normalized if $|T|=1$. We prepare the input state

$\Psi_{\rm
input}=\frac{1}{|T_{0}|^{{n_{bos}}/p}}|T_{0}\rangle^{\otimes{n_{bos}}/p}.$
(42)

Preparing this state takes circuit depth $O(N^{p})$ since we can prepare
${n_{bos}}/p$ copies of the state $\frac{1}{|T_{0}|}|T_{0}\rangle$ in
parallel.

We want to know the expectation value $\langle\Psi_{\rm
input}|E_{PE}|\Psi_{\rm input}\rangle$, but to get oriented, let us estimate
the overlap $\langle\Psi_{\rm input}|\Psi_{\rm sig}\rangle$.

We have

$\langle v_{\rm sig}^{\otimes p}|T_{0}\rangle=\lambda N^{p}+\langle v_{\rm
sig}^{\otimes p}|G\rangle.$ (43)

The probability distribution of $\langle v_{\rm sig}^{\otimes p}|G\rangle$ is
a Gaussian with zero mean and unit variance, so with high probability,
$\langle v_{\rm sig}^{\otimes p}|G\rangle$ is $o(\lambda N^{p})$; indeed, for
any increasing function of $N$ which diverges as $N\rightarrow\infty$, with
high probability it is bounded by that function.

Hence, with high probability,

$\langle v_{\rm
sig}^{\otimes{n_{bos}}}|T_{0}^{\otimes{n_{bos}}/p}\rangle=(1-o(1))\cdot\lambda^{{n_{bos}}/p}N^{n_{bos}}.$
(44)

At the same time,

$\langle
T_{0}^{\otimes{n_{bos}}/p}|T_{0}^{\otimes{n_{bos}}/p}\rangle=|T_{0}|^{2{n_{bos}}/p}.$
(45)

We have $\mathbb{E}[|G|^{2}]=O(N^{p}),$ where the precise constant in the
big-O notation depends on whether we symmetrize $G$ or not and whether we use
complex entries or not. Further, $|G|^{2}$ is a sum of squares of independent
random variables (the entries of $G$). So, by central limit, with high
probability $|G|^{2}$ is bounded by $O(N^{p})$. So, with high probability,
$|T_{0}|^{2{n_{bos}}/p}=O(N^{{n_{bos}}})$.

So, with high probability,

$\frac{\Bigl{|}\langle v_{\rm
sig}^{\otimes{n_{bos}}}|T_{0}^{\otimes{n_{bos}}/p}\rangle\Bigr{|}^{2}}{\langle
v_{\rm sig}^{\otimes{n_{bos}}}|v_{\rm
sig}^{\otimes{n_{bos}}}\rangle\cdot\langle
T_{0}^{\otimes{n_{bos}}/p}|T_{0}^{\otimes{n_{bos}}/p}\rangle}\geq(1-o(1))\cdot\lambda^{2{n_{bos}}/p}.$
(46)

For $\lambda=CN^{-p/4}$, this is $(1-o(1))C^{2{n_{bos}}/p}N^{-{n_{bos}}/2}$.

If $\psi_{\rm target}$ were equal to $\Psi_{\rm sig}=N^{-{n_{bos}}/2}|v_{\rm
sig}^{\otimes{n_{bos}}}\rangle$, then for fixed $N^{-p/4}/\lambda$, Eq. (46)
would give a lower bound to the squared overlap of the initial state with
$\psi_{\rm target}$ which would be quadratically better (in terms of its
scaling with $N$) than the squared overlap for the maximally entangled input
state. So, after applying amplitude amplification, this would give a quartic
improvement over the fastest classical algorithm.

However, $\psi_{\rm target}$ is not equal to $\Psi_{\rm
sig}=N^{-{n_{bos}}/2}|v_{\rm sig}^{\otimes{n_{bos}}}\rangle$ and so this does
not give a lower bound on $\langle\psi_{\rm target}|\Psi_{\rm input}\rangle$.

However, we have:

######  Lemma 8.

Suppose that $E_{0}\geq(1+c)E_{max}$. Then, the projection of $\Psi_{\rm sig}$
onto the eigenspace with eigenvalue $\geq(7/8)E_{0}+(1/8)E_{max}$ is greater
than or equal to $\Omega(1)\cdot c/(1+c)$.

###### Proof.

The expectation value $\langle\Psi_{\rm sig}|H(T_{0})|\Psi_{\rm sig}\rangle$
was estimated in subsections 3.2.1,3.2.2. We have that with high probability,
this expectation value is $\geq(1-o(1))E_{0}$. With high probability, the
largest eigenvalue of $H(T_{0})$ in absolute value is bounded by
$(1+o(1))(E_{0}+E_{max})$; for the even case this is just the triangle
inequality, while for the odd case this uses lemma 3. Hence, by Markov’s
inequality applied to $\lambda_{1}-H(T_{0})$, the projection of $\Psi_{\rm
sig}$ onto the eigenspace with eigenvalue $\geq(7/8)E_{0}+(1/8)E_{max}$ is
greater than or equal to
$(1-o(1))E_{0}-((7/8)E_{0}+(1/8)E_{max})/((1+o(1))(E_{0}+E_{max}))$ which is
$\geq\Omega(1)\cdot c/(1+c)$. ∎

So, $\Psi_{\rm sig}$ has some non-negligible projection onto the desired
eigenspace. This does not however yet give us a lower bound on
$\langle\Psi_{\rm input}|E_{PE}|\Psi_{\rm input}\rangle$: we can expand
$\Psi_{\rm input}$ as a linear combination of $\Psi_{\rm sig}$ and some
orthogonal state but we have not bounded the cross-terms in the expectation
value.

However, we now give heuristic evidence (not a proof) for a lower bound on
$\mathbb{E}[\langle\Psi_{\rm input}|E_{PE}|\Psi_{\rm input}\rangle]$. The main
assumption we will make is that $\lambda_{1}=E_{0}\cdot(1+o(1/\log(N))$; we
expect that that assumption can be proven to hold with high probability. We
consider just the case of even $p$ (we expect that odd $p$ can be handled
similarly).

The main reason that we do not give a full proof is that a lower bound on
$\mathbb{E}[\langle\Psi_{\rm input}|E_{PE}|\Psi_{\rm input}\rangle]$ will only
imply a lower bound on the expected inverse squared runtime, rather than an
upper bound on the expected runtime, which is what we really want. Instead in
the next subsubsection we give a modified algorithm with an upper bound on the
expected runtime. Let us note that we conjecture that algorithm 6 does have a
quartic speedup for the expect runtime with high probability. One might guess
that this could be proven using the bound on expectation value of the inverse
squared runtime and some concentration of measure argument. However, we have
not been able to make this precise.

Let us clarify some terminology to distinguish two meanings of the word
“expected”, corresponding to averages over $G$ or averages over outcomes of a
quantum algorithm, i.e., to the “expected runtime”. From here on, when we
refer to the “runtime” of a phase estimation algorithm, this is a short way of
saying the expected runtime for a given choice of $G$. When we refer to the
“expectation value of the runtime”, we mean the expectation value over $G$ of
this expected runtime. Applying amplitude amplification, the runtime is
bounded by the time for the state preparation and phase estimation multiplied
by the inverse square-root of $\langle\Psi_{\rm input}|E_{PE}|\Psi_{\rm
input}\rangle$. So, lower bounding $\mathbb{E}_{G}[\langle\Psi_{\rm
input}|E_{PE}|\Psi_{\rm input}\rangle]$ will upper bound the expectation value
over $G$ of the inverse squared runtime and we will find a bound on the
expectation value of the inverse squared runtime by

$\Bigl{(}N^{{n_{bos}}/4}{\rm
poly}(N,{n_{bos}},1/(E_{0}-E_{max}),\log(D(N,{n_{bos}})/\epsilon))\Bigr{)}^{-2}\Bigl{(}N^{-p/4}/\lambda\Bigr{)}^{-2{n_{bos}}/p}$

in the case that $E_{0}\geq E_{max}\cdot(1+c)$ for any $c>0$. For fixed
$N^{-p/4}/\lambda$, this gives a further quadratic improvement, in terms of
the scaling of the runtime with $N$, over algorithm 4.

Given the existence of that modified algorithm of subsubsection 5.2.3, we will
just sketch an outline of a possible proof of the lower bound on
$\mathbb{E}[\langle\Psi_{\rm input}|E_{PE}|\Psi_{\rm input}\rangle]$, leaving
several details out. The basic idea is to lower bound this expectation value
by a tensor network, working in some approximation which amounts to ignoring
fluctuations in $\lambda_{1}$, then average the tensor network by a sum of
pairings, and use this sum to lower bound the expectation value. Roughly the
physical idea is that the terms in $\Psi_{\rm input}$ proportional to $G$ will
tend on average to increase the overlap $\langle\Psi_{\rm input}|\psi_{\rm
target}\rangle$, rather than decrease it.

Consider the operator $\lambda_{1}^{-m}H(T_{0})^{m}$ for large $m$. If we take
$m$ sufficiently large compared to ${n_{bos}}\log(N)/\log(E_{0}/E_{max})$,
this will operator will lower bound $E_{PE}$ up to some negligible error. That
is, we take $m$ large enough that $\lambda_{1}^{-m}H(T_{0})^{m}$ is negligibly
small acting on any eigenvector $v_{i}$ with eigenvalue
$\lambda_{i}\leq(7/8)E_{0}+(1/8)E_{max}$, and for
$\lambda_{i}\geq(7/8)E_{0}+(1/8)E_{max}$, Eq. (41) gives a lower bound on
$E_{PE}$ that is equal to $1$ up to some negligible phase estimation error
$\tilde{\epsilon}$ while clearly $\lambda_{1}^{-m}H(T_{0})^{m}\leq 1$. For
fixed $E_{0}/E_{max}$, it suffices to take $m=O({n_{bos}}\log(N))$.

For the range of $m$ that we consider here, by assumption we can ignore the
fluctuations in $\lambda_{1}$, i.e., we approximate

$\displaystyle|\langle\Psi_{\rm input}|E_{PE}|\Psi_{\rm input}\rangle|^{2}$
$\displaystyle\gtrsim$
$\displaystyle(\mathbb{E}[\lambda_{1}])^{-m}\mathbb{E}[\langle\Psi_{\rm
input}|H(T_{0})^{m}|\Psi_{\rm input}\rangle]$ (47) $\displaystyle\approx$
$\displaystyle(\mathbb{E}[\lambda_{1}])^{-m}\frac{1}{\mathbb{E}[\langle
T_{0}^{\otimes{n_{bos}}/p}|T_{0}^{\otimes{n_{bos}}/p}\rangle]}[\mathbb{E}[\langle
T_{0}^{\otimes{n_{bos}}/p}|H(T_{0})^{m}|T_{0}^{\otimes{n_{bos}}}\rangle],$
(48)

where in the second line we further approximate that we can ignore
fluctuations in the norm $\langle
T_{0}^{\otimes{n_{bos}}/p}|T_{0}^{\otimes{n_{bos}}/p}\rangle$ and treat it as
a constant (the proof is a standard large deviation argument on the norm of
the tensor).

The quantity $\langle
T_{0}^{\otimes{n_{bos}}/p}|H(T_{0})^{m}|T_{0}^{\otimes{n_{bos}}/p}\rangle$ can
be evaluated by a sum of tensor networks, using the full Hilbert space, i.e.,
for each of $m$ choices of $i_{1},\ldots,i_{p/2}$ in each of the $m$ factors
of $H(T_{0})$ we have a tensor network). We can then write each tensor network
as a sum of tensor networks, inserting either $\lambda v_{\rm sig}^{\otimes
n}$ or $G$ for each tensor, and we can average these tensor networks over $G$
using the methods of appendix A by summing over pairings. Since every term in
this sum over networks and pairings is positive, if we restrict to some subset
of terms, we get a lower bound. Let us restrict to the terms in which for
$\Psi_{\rm input}$, we choose $\lambda v_{\rm sig}^{\otimes p}$ for every
tensor. For this set of terms, the tensor network computes precisely
$\mathbb{E}[\langle\lambda^{{n_{bos}}/p}v_{\rm
sig}^{\otimes{n_{bos}}}|H(T_{0})^{m}|\lambda^{{n_{bos}}/p}v_{\rm
sig}^{\otimes{n_{bos}}}\rangle]$. This in turn is
$\geq(E_{0}\cdot(1-o(1/m))^{m}|\lambda^{{n_{bos}}/p}v_{\rm
sig}^{\otimes{n_{bos}}}|^{2}$, since $\langle\Psi_{\rm sig}|H(T_{0})|\Psi_{\rm
sig}\rangle\geq E_{0}\cdot(1-o(1/m))$ for $m=O({n_{bos}}\log(N))$. So the
tensor network is lower bounded by $E_{0}^{m}|\lambda^{{n_{bos}}/p}v_{\rm
sig}^{\otimes{n_{bos}}}|^{2}$.

So, with these approximations we have lower bounded
$\mathbb{E}[\langle\Psi_{\rm input}|E_{PE}|\Psi_{\rm input}\rangle]$.

We make some implementation remarks on the algorithm. The algorithm as
described requires measuring whether we are in the symmetric subspace. Note
that the input state $\Psi_{\rm input}$ need not be in the symmetric subspace.
Such a projection can be done for example by phase estimating a Hamiltonian
which is a sum of permutation operators. One can also omit this projection
onto the symmetric subspace since our upper bounds on $H(G)$ holds both in the
full Hilbert space and in the symmetric subspace.

We have considered the case that ${n_{bos}}$ is an integer multiple of $p$. If
${n_{bos}}=kp+l$ for some integers $k,l$ with $0<l<p$, then one can use $l$
ancilla qudits, and prepare an input state which is equal to

$\frac{1}{|T_{0}|^{k}}|T_{0}\rangle^{\otimes k},$

on $kp$ qudits, tensored with a maximally entangled state between the
remaining $l$ qudits and the remaining $l$ ancillas. The idea is that we get
the additional quadratic improvement in overlap on $kp$ of the qudits, and the
remaining $l$ ancilla only cost $1/{\rm poly}(N)$ overlap since $l=O(1)$.

####  5.2.3 Chosen Input State: Modified Version

We now modify the algorithm 5 (and its amplified version) to obtain an
algorithm for which we can prove the quadratic improvement over algorithm 4
without any assumption. Consider given $T_{0}$. Let $\Delta$ be a $p$-th order
tensor, chosen from the same distribution as $G$. Consider the tensor

$\displaystyle T_{0}^{\prime}$ $\displaystyle=$ $\displaystyle T_{0}+x\Delta$
$\displaystyle\equiv$ $\displaystyle\lambda v_{\rm sig}^{\otimes
p}+G^{\prime}$

for some real scalar $x$, where the tensor $G^{\prime}\equiv G+x\Delta$ has
Gaussian entries with variance of the entries equal to $1+x^{2}$. We will
assume $x=O(1)$; indeed later we will choose $x=o(1)$. Let us write $\Psi_{\rm
input}(T)\equiv|T|^{-{n_{bos}}/p}|T^{\otimes{n_{bos}}/p}\rangle$ and
$E_{PE}(T)$ to denote the phase estimation operator $E_{PE}$ for Hamiltonian
$H(T)$.

We have

$\displaystyle G$ $\displaystyle=$
$\displaystyle\frac{1}{1+x^{2}}(G+x\Delta)+\frac{x}{1+x^{2}}(xG-\Delta)$
$\displaystyle=$
$\displaystyle\frac{1}{1+x^{2}}G^{\prime}+\frac{x}{\sqrt{1+x^{2}}}\frac{(xG-\Delta)}{\sqrt{1+x^{2}}}$
$\displaystyle=$
$\displaystyle\frac{1}{1+x^{2}}G^{\prime}+\frac{x}{\sqrt{1+x^{2}}}\delta,$

where $\delta=(1+x^{2})^{-1/2}(xG-\Delta)$. The two random variables
$G^{\prime}$ and $\delta$ have vanishing covariance, so Eq. (5.2.3) expresses
$G$ as a scalar multiple of $G^{\prime}$ plus an additional Gaussian random
variable $\delta$ which is independent of $G^{\prime}$. The variable $\delta$
also has variance $1$.

Given $G$ and $G^{\prime}$, let

$G(y)=yG^{\prime}+(1-y)G=(y+\frac{1-y}{1+x^{2}})G^{\prime}+\frac{x(1-y)}{\sqrt{1+x^{2}}}\delta.$
(51)

so that $G(y)$ linearly interpolates between $G$ and $G^{\prime}$.

The idea behind the algorithm is to take a given $G$ as input, randomly
perturb to produce $G^{\prime}$, and then consider several input states
$\Psi_{\rm input}(\lambda v_{\rm sig}^{\otimes p}+G(y))$ with different
choices of $y\in[0,1]$.

Let us first recall a property of normalization. We have $\Psi_{\rm
input}=\frac{1}{|T_{0}|^{{n_{bos}}/p}}|T_{0}\rangle^{\otimes{n_{bos}}/p}.$ Let
us write

$Z=|T_{0}|^{2{n_{bos}}/p},$

so $\Psi_{\rm input}=Z^{-1/2}|T_{0}\rangle^{\otimes{n_{bos}}/p}.$ As shown
before, with high probability
$Z^{1/2}=|T_{0}|^{{n_{bos}}/p}=O(N^{{n_{bos}}/2})$. Further, with high
probability the fluctuations of $Z^{1/2}$ are $o(1)$ compared to its
expectation value. So, from here on we will treat this normalization factor
$Z$ as a constant, i.e., of course the normalization depends on $N,{n_{bos}}$
but we will ignore its dependence on $G$. We emphasize that we are not making
any additional assumption here as with high probability the fluctuations are
asymptotically negligible; we are simply choosing not to write the
normalization explicitly. (Remark: indeed, all we really need is an upper
bound on $|T_{0}|^{{n_{bos}}/p}$ that holds with high probability, since the
normalization constant $|T_{0}|^{{n_{bos}}/p}$ always appears in the
denominator.) Further, we will, without remarking on it further, treat other
normalization factors such as $|\Psi_{\rm input}(\lambda v_{\rm sig}^{\otimes
p}+G(y))|$ as constants, and we will introduce other notation for those
constants. Indeed, because we treat the normalization factors such as $Z$ as
constants, we will mostly work with un-normalized states which simplifies some
of the calculations.

In an abuse of notation, let us define $\Psi_{\rm input}(y)=\Psi_{\rm
input}(\lambda v_{\rm sig}^{\otimes p}+G(y))$, and write $Z(y)=|\lambda v_{\rm
sig}^{\otimes p}+G(y)|^{2{n_{bos}}/p}$. Let us write $T_{0}(y)=\lambda v_{\rm
sig}^{\otimes p}+G(y)$.

Let $\Psi(y)$ denote the un-normalized state
$|T_{0}(y)\rangle^{\otimes{n_{bos}}/p}$ so that $\Psi_{\rm
input}(y)=Z(y)^{-1/2}\Psi(y)$. Let us expand $\Psi(y)$ as a series in $\delta$
and define $\Psi^{0}(y)$ to denote the zeroth order term in $\delta$.

As a warmup, let us consider
$\langle\Psi^{0}(y)|E_{PE}(T_{0}^{\prime})|\Psi^{0}(y)\rangle$. We consider
the higher order terms in $\delta$ later.

We expand the state $\Psi^{0}(y)$ as a series in $(y+\frac{1-y}{1+x^{2}})$.
Doing this means that we express
$\langle\Psi^{0}(y)|E_{PE}(T_{0}^{\prime})|\Psi^{0}(y)\rangle$ as a polynomial
of degree $2{n_{bos}}$ which we write as $\sum_{i\geq
0}a_{i}(y+\frac{1-y}{1+x^{2}})^{i}.$ The zero-th order term $a_{0}$ is simply
equal to $\langle(\lambda v_{\rm sig}^{\otimes
p})^{\otimes{n_{bos}}/p}|E_{PE}(T_{0}^{\prime})|(\lambda v_{\rm sig}^{\otimes
p})^{\otimes{n_{bos}}/p}\rangle$, which is lower bounded in lemma 8. Hence,

$a_{0}\geq\lambda^{2{n_{bos}}/p}N^{{n_{bos}}}\Omega(1)\cdot c/(1+c)$ (52)

for $E_{0}\geq(1+c)E_{max}$.

Now we use a result about polynomials:

######  Lemma 9.

Let $p(z)$ be a polynomial of degree $2{n_{bos}}$. Let $[a,b]$ be an interval
with $0\leq a<b$. Then,

${\rm
max}_{z\in[a,b]}|p(z)|\geq\Bigl{|}\frac{a+b}{b-a}\Bigr{|}^{-2{n_{bos}}}\exp(-O({n_{bos}}))|p(0)|.$
(53)

###### Proof.

We minimize ${\rm max}_{z\in[a,b]}|p(z)|$ over polynomials of degree
$2{n_{bos}}$ with given value $p(0)$. Applying an affine transformation
$z\rightarrow(2/(b-a))(z-(a+b)/2)$, this is equivalent to minimizing ${\rm
max}_{z\in[-1,1]}|p(z)|$ over polynomials of degree $2$ with given value
$p(z_{0})$ for $z_{0}=(a+b)/(a-b)$. We claim that this is minimized by

$\frac{p(z_{0})}{T_{2{n_{bos}}}(z_{0})}T_{2{n_{bos}}}(z),$

where $T_{2{n_{bos}}}$ is a Chebyshev polynomial. Proof of claim: suppose some
other polynomial $q(z)$ has a smaller maximum absolute value on $[-1,+1]$ with
$q(z_{0})=p(z_{0})$. Then the polynomial $p(z)-q(z)$ has a zero at $z_{0}$ but
also has at least $2{n_{bos}}$ zeros on the interval $[-1,+1]$; this follows
from the intermediate value theorem because $T_{2{n_{bos}}}$ has
$2{n_{bos}}+1$ extreme points on the interval which alternate signs. This
gives a contradiction since $p(z)-q(z)$ is degree at most $2{n_{bos}}$.

So, ${\rm max}_{z\in[a,b]}|p(z)|\geq|p(z_{0})|/|T_{2{n_{bos}}}(z_{0})|$. If
$0\not\in[a,b]$ then $|z_{0}|>1$. We can bound $T_{2{n_{bos}}}(z_{0})$ for
$|z_{0}|>1$ by $|z_{0}|^{2{n_{bos}}}$ times the sum of absolute values of
coefficients of $T_{2{n_{bos}}}$. This sum of coefficients is bounded by
$\exp(O({n_{bos}}))$. ∎

We apply this lemma to lower bound ${\rm max}_{y\in[0,1]}|\sum_{i\geq
0}a_{i}(y+\frac{1-y}{1+x^{2}})^{i}|$. For $z=y+\frac{1-y}{1+x^{2}}$, with
$b=1$ and $a=(1+x^{2})^{-1}$, this is ${\rm max}_{z\in[a,b]}|p(z)|$ with
$p(z)=\sum_{i}a_{i}z^{i}$. So, using that $(a+b)/(b-a)\leq(2/x^{2})$, we have

${\rm max}_{y\in[0,1]}|\sum_{i\geq 0}a_{i}(y+\frac{1-y}{1+x^{2}})^{i}|\geq
x^{4{n_{bos}}}\exp(-O({n_{bos}}))|a_{0}|.$ (54)

However, for purposes of an algorithm, we need to consider not just the
maximum of this polynomial $p(z)$ over an interval, but also consider whether
we can come close to this maximum by sampling it at some small number of
discrete points. Fortunately we have the following lemma:

######  Lemma 10.

Let $p(z)$ be a polynomial of degree $2{n_{bos}}$. Let $[a,b]$ be an interval
which does not contain $0$. Then, with probability at least $1/2$, if one
selects a random point $z$ in the interval $[a,b]$ from the uniform
distribution, we have
$|p(z)|\geq\Bigl{|}\frac{a+b}{b-a}\Bigr{|}^{-2{n_{bos}}}\exp(-O({n_{bos}}))|p(0)|$.

###### Proof.

Apply an affine transformation $z\rightarrow(2/(b-a))(z-(a+b)/2)$, which maps
$a$ to $-1$, $b$ to $+1$ and $0$ to $(a+b)/(a-b)$. Write
$p(z)=\prod_{i}(z-z_{i})$ where $z_{i}$ are zeros of the polynomial; we
multiply by $p(z)$ by a normalization constant so that the highest order
coefficient is equal to $1$. Then $\log(|p(z)|)=\sum_{i}{\rm
Re}(\log(z-z_{i}))$. Let $A$ denote the average of this logarithm over the
interval $[-1,+1]$; by calculus this is

$A=\frac{1}{2}\sum_{i}{\rm
Re}\Bigl{(}(z-z_{i})\log(z-z_{i})-(z-z_{i})\Bigr{)}\Bigl{|}_{z=-1}^{z=+1}.$

We claim that

$\log(p((a+b)/(a-b)))-A\leq
2{n_{bos}}\cdot\Bigl{(}\log((a+b)/(a-b))+O(1)\Bigr{)}.$ (55)

To show this let $T_{i}$ denote a given term in the sum for $A$, i.e.
$T_{i}=(1/2){\rm Re}((z-z_{i})\log(z-z_{i})-(z-z_{i}))\Bigl{|}_{z=-1}^{z=+1}$.
Let $D_{i}\equiv{\rm Re}(\log((a+b)/(a-b)-z_{i}))-T_{i}$. By considering
various cases, we will show that $D_{i}\leq\log((a+b)/(a-b))+O(1)$ which
implies Eq. (55). Proof of claim: for $|z_{i}|$ less than or equal to some
fixed constant (for example, $|z_{i}|\leq 10$), $T_{i}$ is lower bounded by
some absolute constant, so $D_{i}\leq\log((a+b)/(a-b))+O(1)$. For $|z_{i}|$
larger than this fixed constant (for example, $|z_{i}|>10$), $T_{i}$ is at
least $\log(|z_{i}|)$ minus some other absolute constant, so $D_{i}\leq{\rm
Re}(\log((a+b)/(a-b)-z_{i})-\log(z_{i}))+O(1)={\rm
Re}(\log((a+b)/(a-b)))+O(1)$.

Further, for each $i$, we see that ${\rm max}_{z\in[-1,+1]}{\rm
Re}(\log(z-z_{i}))$ is upper bounded by $T_{i}+O(1)$. Hence, ${\rm
max}_{z\in[-1,+1]}\log(|p(z)|)$ is upper bounded by $A+O({n_{bos}})$. Hence,
with probability at least $1/2$, for $z$ chosen uniformly in $[-1,+1]$, we
have that $\log(|p(z)|)\geq A-O({n_{bos}})$.

Hence, with probability at least $1/2$, we have that $\log(|p(z)|)\geq
A-2{n_{bos}}\cdot(\log((a+b)/(a-b))-O({n_{bos}})$. ∎

So, noting that for $z=y+\frac{1-y}{1+x^{2}}$, a uniform choice of $y$ on
$[0,1]$ is the same as a uniform choice of $z$ on $[1/(1+x^{2}),1]$, we have

######  Lemma 11.

For $y$ chosen randomly from the uniform distribution on $[0,1]$, the quantity
$\langle\Psi^{0}(y)|E_{PE}(T_{0}^{\prime})|\Psi^{0}(y)\rangle$ is greater than
or equal to $\langle(\lambda v_{\rm sig}^{\otimes
p})^{\otimes{n_{bos}}/p}|E_{PE}(T_{0}^{\prime})|(\lambda v_{\rm sig}^{\otimes
p})^{\otimes{n_{bos}}/p}\rangle\cdot\exp(-O({n_{bos}}))x^{4{n_{bos}}}$.

Using Eq. (52),

$\langle\Psi^{0}(y)|E_{PE}(T_{0}^{\prime})|\Psi^{0}(y)\rangle\geq\lambda^{2{n_{bos}}/p}N^{{n_{bos}}}\exp(-O({n_{bos}}))x^{4{n_{bos}}}\cdot
c/(1+c).$ (56)

We later will worry about optimizing over $x$. Increasing $x$ will increase
the right-hand side of Eq. (56), but it will also change the normalization
constant $Z(y)$ that we must multiply by to correctly normalize the input
state and also will worsen the threshold for recovery (since $E_{max}$ will
increase). For now, let us note for orientation that if we choose, for
example, $x=0.01$, then the tensor $T_{0}^{\prime}$ is chosen from the same
distribution as $T_{0}$ up to a multiplication by a factor of
$(1+0.01^{2})^{1/2}$ which is very close to $1$. This leads to only a very
small reduction in the threshold $\lambda$ at which recovery is possible. At
the same time, the change in normalization $(1+0.01^{2})^{{n_{bos}}}$ is
asymptotically negligible (for fixed ${n_{bos}}$) compared to the improvement
which is polynomial in $N^{n_{bos}}$.

We now consider the additional terms depending on $\delta$, i.e., due to the
difference $\Psi(y)-\Psi^{0}(y)$. At first glance, we have not gained much by
the above trick introducing $T^{\prime}$ since we still have this Gaussian
random variable $\delta$ to consider. However, we have the advantage now that
in $\langle\Psi(y)|E_{PE}(T_{0}^{\prime})|\Psi(y)\rangle$, the operator
$E_{PE}(T_{0}^{\prime})$ does not depend on $\delta$ and the states in the bra
and ket depend on $\delta$ only polynomially (treating the overall
normalization as a constant).

Let us outline our approach. The main worry that we have is that this quantity
$\langle\Psi(y)|E_{PE}(T_{0}^{\prime})|\Psi(y)\rangle$ will have a probability
distribution that is peaked near some “small” value, i.e., much less than the
value of the overlap at $\delta=0$, rather than some “large” value, i.e.,
roughly comparable to the value of the overlap at $\delta=0$. To deal with
this, we will use a trick of adding additional noise to make the expectation
value of the overlap large, i.e., roughly comparable to the value of the
overlap at $\delta=0$. Now we still need to worry about whether the
probability function might have some large probability of having a small
value. So, we will then appeal we appeal to a theorem of Carbery-Wright on
“anti-concentration” of polynomials of Gaussian random variables, bounding the
probability that the probability lies in some small interval. This theorem
gives us a useful bound unless the variance of the polynomial is small;
however, in that case we can show that the polynomial is likely to be close to
its expectation value.

The trick of adding additional noise is as follows. We perturb the input state
by adding additional Gaussian random noise. Let $\Psi(y,x^{\prime})$ denote
the unnormalized state

$|\lambda v_{\rm sig}^{\otimes
p}+(y+\frac{1-y}{1+x^{2}})G^{\prime}+\frac{x(1-y)}{\sqrt{1+x^{2}}}(\delta+\delta^{\prime})\rangle^{\otimes{n_{bos}}/p},$

where $\delta^{\prime}$ is an additional tensor chosen from a Gaussian
distribution with some standard deviation $x^{\prime}$. So, the tensor
$\delta+\delta^{\prime}$ is sampled from some distribution with variance
$1+x^{\prime 2}$. Let $Z(y,x^{\prime})=|\Psi(y,x^{\prime})|^{2}$. Let
$\Psi_{\rm input}(y,x^{\prime})=Z(y,x^{\prime})^{-1/2}\Psi(y,x^{\prime})$ be
the “perturbed input state”.

Then, we consider the expectation value over $\delta^{\prime}$ of
$\langle\Psi(y,x^{\prime})|E_{PE}(T_{0}^{\prime})|\Psi(y,x^{\prime})\rangle$;
we consider this expectation value as a series in the variable $1+x^{\prime
2}$. Using the same treatment as above of this expectation value as a
polynomial, we find that if $x^{\prime 2}$ is chosen uniformly from the
interval $[0,x^{2}]$ then, with probability at least $1/2$, the quantum
expectation value
$\langle\Psi(y,x^{\prime})|E_{PE}(T_{0}^{\prime})|\Psi(y,x^{\prime})\rangle$
is at least equal to the quantum expectation value with $\delta^{\prime}=0$
multiplied by $\exp(-O({n_{bos}}))x^{8{n_{bos}}}$.

Hence,

######  Lemma 12.

With probability at least $1/4$, for uniform random choices of $y$ and
$x^{\prime 2}$ on $[0,1]$ and $[0,x^{2}]$, the expectation value over
$\delta,\delta^{\prime}$ of
$\langle\Psi(y,x^{\prime})|E_{PE}(T_{0}^{\prime})|\Psi(y,x^{\prime})\rangle$
is at least $\langle(\lambda v_{\rm sig}^{\otimes
p})^{\otimes{n_{bos}}/p}|E_{PE}(T_{0}^{\prime})|(\lambda v_{\rm sig}^{\otimes
p})^{\otimes{n_{bos}}/p}\rangle\cdot\exp(-O({n_{bos}}))x^{8{n_{bos}}}$.

We emphasize that throughout we are treating $G^{\prime}$ as fixed and
considering $\delta,\delta^{\prime}$ as random. Now we use anti-concentration
in the Gaussian random variables $\delta,\delta^{\prime}$. The following lemma
is a consequence of the Carbery-Wright theorem [21]

######  Lemma 13.

Let $p(z_{1},\ldots,z_{n})$ be a polynomial of degree at most $k$ in $n$
independent normally distributed random variables. Then, for any
$t\in{\mathbb{R}}$ and any $\delta>0$ with $|t-\mathbb{E}[p]|>\delta$,

$\displaystyle{\rm Pr}[|p(z_{1},\ldots,z_{n})-t|\leq\delta]$
$\displaystyle\leq$
$\displaystyle\Bigl{(}\frac{\delta}{|t-\mathbb{E}[p]|-\delta}\Bigr{)}^{2/(2k+1)}O(k)^{2k/(2k+1)}$
$\displaystyle\leq$
$\displaystyle\Bigl{(}\frac{\delta}{|t-\mathbb{E}[p]|-\delta}\Bigr{)}^{2/(2k+1)}O(k),$

where the big-O notation here hides a universal constant.

As a corollary, choosing $t=0$,

${\rm
Pr}[|p|\leq\delta|\mathbb{E}[p]|]\leq\Bigl{(}\frac{\delta}{1-\delta}\Bigr{)}^{2/(2k+1)}O(k)=O(\delta^{2/(2k+1)})O(k).$
(58)

###### Proof.

Let ${\rm Var}(p)$ denote the variance of $p(\cdot)$. As a trivial bound,

${\rm Pr}[|p(z_{1},\ldots,z_{n})-t|\leq\delta]\leq{\rm
Var}(p)/\Bigl{(}|t-\mathbb{E}[p]|-\delta\Bigr{)}^{2}.$ (59)

By Carbery-Wright,

${\rm Pr}[|p(z_{1},\ldots,z_{n})-t|\leq\delta]\leq O(k)\cdot(\delta/\sqrt{{\rm
Var}(p)})^{1/k}.$ (60)

Maximizing the bound from Eqs. (59,60) over ${\rm Var}(p)$, Eq. (13) follows.
∎

Hence, from Eq. (58),

######  Lemma 14.

Consider a choice of $y,x^{\prime}$ such that the expectation value over
$\delta,\delta^{\prime}$ of the quantum expectation value
$\langle\Psi(y,x^{\prime})|E_{PE}(T_{0}^{\prime})|\Psi(y,x^{\prime})\rangle$
is at least $\langle(\lambda v_{\rm sig}^{\otimes
p})^{\otimes{n_{bos}}/p}|E_{PE}(T_{0}^{\prime})|(\lambda v_{\rm sig}^{\otimes
p})^{\otimes{n_{bos}}/p}\rangle\cdot\exp(-O({n_{bos}}))x^{8{n_{bos}}}$. Then,
for random choice of $\delta,\delta^{\prime}$, the quantum expectation value
$\langle\Psi(y,x^{\prime})|E_{PE}(T_{0}^{\prime})|\Psi(y,x^{\prime})\rangle$
is at least

$\langle(\lambda v_{\rm sig}^{\otimes
p})^{\otimes{n_{bos}}/p}|E_{PE}(T_{0}^{\prime})|(\lambda v_{\rm sig}^{\otimes
p})^{\otimes{n_{bos}}/p}\rangle\cdot\exp(-O({n_{bos}}))x^{8{n_{bos}}}\delta,$

with probability at least $1-O({n_{bos}})O(\delta)^{2/(2{n_{bos}}+1)}$.

Lemma 14 considers unnormalized input states. Dividing by normalization
constant $Z(y,x^{\prime})$ we can get a lower bound on the expectation value
for the normalized input state $\Psi_{\rm input}(y,x^{\prime})$ by (for
$E_{0}\geq(1+c)E_{max}$):

$\displaystyle\frac{1}{Z(y,x^{\prime})}\langle(\lambda v_{\rm sig}^{\otimes
p})^{\otimes{n_{bos}}/p}|E_{PE}(T_{0}^{\prime})|(\lambda v_{\rm sig}^{\otimes
p})^{\otimes{n_{bos}}/p}\rangle\cdot\exp(-O({n_{bos}}))x^{8{n_{bos}}}\delta$
$\displaystyle=$
$\displaystyle\lambda^{2{n_{bos}}/p}\exp(-O({n_{bos}}))x^{8{n_{bos}}}\delta\cdot\Omega(1)$
$\displaystyle=$
$\displaystyle\Bigl{(}\frac{\lambda}{N^{-p/4}}\Bigr{)}^{2{n_{bos}}/p}N^{-{n_{bos}}/2}x^{8{n_{bos}}}\delta\cdot\Omega(1)$
$\displaystyle=$
$\displaystyle\Bigl{(}\frac{\lambda}{N^{-p/4}}\Bigr{)}^{2{n_{bos}}/p}D(N,{n_{bos}})^{-1/2}x^{8{n_{bos}}}\delta\cdot\Omega(1).$

Finally, let us pick $x=1/\log(N)$. Then, we consider the following algorithm
7.

Algorithm 7 Quantum Algorithm (modified improved input state, unamplified
version)

  * 1.

Choose random $y,x^{\prime}$. Sample $T_{0}^{\prime},\delta^{\prime}$
randomly. Prepare input state $\Psi_{\rm input}(y,x^{\prime})$.

  * 2.

If the initial state is not in the symmetric subspace, report “failure”. If
the state is in the symmetric subspace, apply phase estimation using
$H(T_{0})$ Let $\psi$ be the resulting state. If the resulting eigenvalue is
larger than $(E_{0}+E_{cut})/2$, report “success”. Otherwise, report
“failure”.

  * 3.

If success is reported, measure and return
$\langle\psi|a^{\dagger}_{\mu}a_{\nu}|\psi\rangle.$

We apply amplitude amplification to algorithm 7. We apply amplitude
amplification under the assumption that indeed
$\langle\Psi(y,x^{\prime})|E_{PE}(T_{0}^{\prime})|\Psi(y,x^{\prime})\rangle$
is at least $\langle(\lambda v_{\rm sig}^{\otimes
p})^{\otimes{n_{bos}}/p}|E_{PE}(T_{0}^{\prime})|(\lambda v_{\rm sig}^{\otimes
p})^{\otimes{n_{bos}}/p}\rangle\cdot\exp(-O({n_{bos}}))x^{4{n_{bos}}}\delta.$
From the theorems above, with high probability in $\delta^{\prime}$, this
happens at least $1/4$ of the time. If the assumption does not hold, we re-
sample $y,x^{\prime}$ and try again. Then, we find that for $E_{0}\geq