E_{max}\cdot(1+O(1/\log(N))$, with probability at least
$1-O({n_{bos}})\delta^{2/(2{n_{bos}}+1)}$, the expected runtime of the
algorithm is at most

${\rm
poly}(N,{n_{bos}},1/(E_{0}-E_{max},\log(D(N,{n_{bos}})/\epsilon))\exp(O({n_{bos}}))\log(N)^{4{n_{bos}}}(\frac{N^{-p/4}}{\lambda})^{{n_{bos}}/p}D(N,{n_{bos}})^{1/4}\delta^{-1/2}.$

Picking $\delta$ to be a slowly decaying function which is $o(1)$, such as
$\delta=1/\log(N)$, we find that

######  Theorem 6.

Let Assumption 1 hold. For $E_{0}\geq E_{max}\cdot(1+c)$, for any $c>0$, with
high probability, the expected runtime of the algorithm is at most

${\rm
poly}(N,{n_{bos}},1/(E_{0}-E_{max},\log(D(N,{n_{bos}})/\epsilon))\exp(O({n_{bos}}))\log(N)^{4{n_{bos}}}(\frac{N^{-p/4}}{\lambda})^{{n_{bos}}/p}D(N,{n_{bos}})^{1/4}.$

Remark: the reader may wonder why we picked $x=1/\log(N)$ rather than a more
slowly decaying function. Indeed, any choice of $x$ which is $o(1)$ would
suffice, and this would replace the factor $\log(N)^{4{n_{bos}}}$ with some
other slowly decaying function.

###  5.3 Further Improvements

In this section, we sketch one improvement that might be practically useful
even if it is not theoretically interesting. For that reason, we do not give
any formal proofs in this section. As we have explained the algorithm above,
once the algorithm succeeds in projecting onto the leading eigenvector, it
measures some element of the single particle density matrix. It is necessary
to repeat this operation, including both preparing the leading eigenvector and
measuring, a total of ${\rm poly}(N,1/\log(\epsilon))$ times to measure the
single particle density matrix to error $\epsilon$, and, as explained, the
algorithm needs to try again to prepare the leading eigenvector each time. The
multiplicative overhead ${\rm poly}(N,1/\log(\epsilon))$ can be turned into an
additive overhead using the non-destructive measurement technique introduced
in Ref. [22]. To implement this, we work in the full Hilbert space (the
leading eigenvector preparation can be implemented in the symetric subspace,
but after preparing the leading eigenvector we work in the full Hilbert
space). Then, we implement a two-outcome projective measurement on a single
qudit (this measurement takes us outside the symmetric subspace though it does
leave us in a subspace which is symmetric under permutation of a subset of
${n_{bos}}-1$ of the qudits), such as $|\mu\rangle_{1}\langle\mu|$ for some
$\mu$. Using the non-destructive measurement techniques of Ref. [22], we can
alternate this two-outcome projective measurement with a projector onto the
leading eigenvector (implemented approximately using phase estimation as
above) to recover the leading eigenvector. See Ref. [22] for more details.

##  6 Networks

In this paper, we have used the language of tensor networks to perform some of
the estimates. In this section, we explore the relation to tensor networks
further. Let us consider the simplest task of detection so that the spectral
algorithms need only to estimate the largest eigenvalue of $H(T_{0})$. In this
case, it suffices to evaluate ${\rm tr}(H(T_{0})^{m})$ to sufficiently large
power $m$. This trace can be written as a sum of tensor networks: the sum is
present because $H(T_{0})$ is a sum of terms. If we restrict to the minimum
${n_{bos}}$ (i.e. ${n_{bos}}=p/2$ for even $p$ and ${n_{bos}}=p-1$ for odd
$p$) then no such sum is present.

Consider a given tensor network $N(T)$ in the sum with $N_{T}$ tensors in the
network. We write the network $N(T)$ to indicate that the graph and ordering
of vertices is given, but that the tensor corresponding to each vertex will be
determined by some tensor $T$ (for some vertices it will be $T$ and for some
it will be $\overline{T}$). We have $N_{T}=m$ for even $p$ or $2m$ for odd
$p$. Thus, the network has $N_{E}=pN_{T}/2$ edges. The tensor network
evaluates to some scalar that we write ${\rm Val}(N(T))$. To get oriented,
suppose that we insert $\lambda v_{\rm sig}^{\otimes p}$ for every tensor in
the network. In this case, the scalar ${\rm Val}(N(\lambda v_{\rm
sig}^{\otimes p}))$ is equal to $\lambda^{N_{T}}N^{N_{E}}$. Suppose instead
that we insert $G$ or $\overline{G}$ for every tensor in the network. It is
possible that $\mathbb{E}[{\rm Val}(N)]$ is zero and it is possible that it is
nonzero. However, let us compute $\mathbb{E}[|{\rm Val}(N(G))|^{2}]$. This
expectation value can be evaluated using the method of appendix A. Note that
${\rm Val}(N(G))|^{2}$ is equal to the evaluation of some tensor network with
$2N_{T}$ tensors; this network is given by two copies of $N$, with all the
tensors complex conjugated in one of the two copies. We refer to the two
copies as the first and second copy. We evaluate this network with $2N_{T}$
tensors by summing over pairings. Consider just the pairing in which we pair
each tensor in the first copy with the corresponding tensor in the second
copy. It is easy to see that this pairing contributes $N^{N_{E}}$ to the
expectation value. Hence, since all pairings contribute positively to the
expectation value, we have that $\mathbb{E}[|{\rm Val}(N)|^{2}]^{1/2}\geq
N^{N_{E}/2}$.

So, we see that $\mathbb{E}[|{\rm Val}(N(G))|^{2}]^{1/2}\gg{\rm Val}(N(\lambda
v_{\rm sig}^{\otimes p}))$ for $N^{N_{E}/2}\gg\lambda^{N_{T}}N^{N_{E}}$, i.e.,
for $\lambda\ll N^{-p/4}$. However, this does not yet imply that
$\mathbb{E}[|{\rm Val}(N(G))|^{2}]^{1/2}\gg\mathbb{E}[{\rm Val}(N(T_{0}))]$,
because, in particular, $\mathbb{E}[{\rm Val}(N(T_{0}))]$ may be larger than
${\rm Val}(N(\lambda v_{\rm sig}^{\otimes p}))$. However, for any network $N$
suppose we insert $\lambda v_{\rm sig}^{\otimes p}$ on some set of tensors $S$
in the network and insert $G$ for the remaining tensors in the network. In
this case, the expectation value is some number given by summing over pairings
of the $N_{T}-q$ insertions of $G$. For a given pairing $P$, consider the
following pairing $Q$ on the tensor network which computes $\mathbb{E}[|{\rm
Val}(N(G))|^{2}]$. For each pair of tensors $T_{1},T_{2}$ in $P$, we pair the
first copy of $T_{1}$ with the first copy of $T_{2}$ and we pair the second
copy of $T_{1}$ with the second copy of $T_{2}$. On the other hand, if a
tensor $T$ is not paired in $P$ (so $T$ is in $S$), we pair the first copy of
$T$ with the second copy of $T$. Using this pairing, we again find, for
$\lambda\ll N^{-p/4}$, for the given choices of inserting $\lambda v_{\rm
sig}^{\otimes p}$ and $G$ into tensors in the network for $\mathbb{E}[{\rm
Val}(N(T_{0}))]$, that the contribution of pairing $Q$ to $\mathbb{E}[|{\rm
Val}(N(G))|^{2}]$ is much larger than the square of the contribution of
pairing $P$ to $\mathbb{E}[{\rm Val}(N(T_{0}))]$.

So, in this regime, the rms fluctuations in the expectation value are large
compared to the expectation value itself, i.e., it seems like detection is not
possible. So, to what extent does this result on tensor networks imply that
spectral algorithms cannot perform inference for $\lambda<N^{-p/4}$? Of
course, such inference is possible by taking ${n_{bos}}$ large so certainly no
general impossibility statement can be made. One way in which one can evade
this result is that the result is for a specific tensor network, while the
algorithms compute a sum over tensor networks. Hence, it is possible that by
taking a large number of networks in the sum, the rms fluctuations in the sum
over networks are small compared to the expectation value of the sum of tensor
networks. This seems to be what happens for the specific algorithms here by
taking ${n_{bos}}$ large giving a large number of networks in the sum. This
connection to tensor networks may be worth exploring further.

Indeed, the way that we were led to develop the spectral algorithm in this
paper is by thinking about different tensor network contractions that might
encode useful information for solving inference problems. After all, a large
class of algorithms can be expressed as evaluating in some way a tensor
network, and there are many circumstances in which a quantum computer may have
some advantage for evaluating a tensor network (given a large number of
caveats, for example acting on some vector encoded as a quantum state by a
tensor encoded as a quantum operator may be a non-unitary operation which may
require some post-selection to implement). Then, a natural class of tensor
networks (in a sense, the “most symmetric”) are those describing the
Hamiltonians considered here. By considering other tensor networks, it may be
possible to apply these methods to other inference problems.

##  7 Discussion

The main new results of this paper are a sequence of spectral algorithms for
odd $p$, and a quantum algorithm to compute the single particle density matrix
of some vector in an eigenspace with large eigenvalue. The quantum algorithm
could be applied to other matrices other than $H(T_{0})$ considered here such
as $Y$ from Ref. [1]; it is not certain though if the speedup using an input
state chosen based on $T_{0}$ would apply to such other matrices.

In this paper, we have related the spectrum for odd $p$ to that for even
$q=2(p-1)$, up to an overall scale factor $\sqrt{N}$, at least in the case of
Gaussian distributed entries of $G$. This enables us to prove correctness of
the algorithm for odd $p$. Also, for ${n_{bos}}=p-1$, the spectrum for
$q=2(p-1)$ is given by random matrix theory, removing a logarithmic factor
present for odd $p$ in some previous work. See corollary 2. We have not
considered the case of entries of $T_{0}$ being chosen from a discrete
distribution, however, we believe that a similar reduction from odd $p$ to
even $q=2(p-1)$ could be proven in that case.

We have used entries of $G$ chosen from a complex Gaussian distribution in the
case of odd $p$ to simplify the proof of the bounds on the spectrum for odd
$p$. We conjecture however that a modification discussed above to cancel
certain “pairing terms” would work for real entries of $G$.

## References

  * [1] Alexander S Wein, Ahmed El Alaoui, and Cristopher Moore.  The kikuchi hierarchy and tensor pca.  2019\.  arXiv:1904.03858.
  * [2] Andrea Montanari and Emile Richard.  A statistical model for tensor pca.  In Advances in Neural Information Processing Systems, pages 2897–2905, 2014.
  * [3] Thibault Lesieur, Leo Miolane, Marc Lelarge, Florent Krzakala, and Lenka Zdeborova.  Statistical and computational phase transitions in spiked tensor estimation.  In 2017 IEEE International Symposium on Information Theory (ISIT). IEEE, jun 2017.  doi:10.1109/isit.2017.8006580.
  * [4] Samuel B Hopkins, Jonathan Shi, and David Steurer.  Tensor principal component analysis via sum-of-square proofs.  In Conference on Learning Theory, pages 956–1006, 2015.
  * [5] Samuel B. Hopkins, Tselil Schramm, Jonathan Shi, and David Steurer.  Fast spectral algorithms from sum-of-squares proofs: tensor decomposition and planted sparse vectors.  In Proceedings of the 48th Annual ACM SIGACT Symposium on Theory of Computing - STOC 2016. ACM Press, 2016.  doi:10.1145/2897518.2897529.
  * [6] Vijay V. S. P. Bhattiprolu, Mrinalkanti Ghosh, Euiwoong Lee, Venkatesan Guruswami, and Madhur Tulsiani.  Multiplicative approximations for polynomial optimization over the unit sphere.  CoRR, abs/1611.05998, 2016.  URL: http://arxiv.org/abs/1611.05998, arXiv:1611.05998.
  * [7] Vijay Bhattiprolu, Mrinalkanti Ghosh, Venkatesan Guruswami, Euiwoong Lee, and Madhur Tulsiani.  Weak decoupling, polynomial folds and approximate optimization over the sphere.  In 2017 IEEE 58th Annual Symposium on Foundations of Computer Science (FOCS). IEEE, oct 2017.  doi:10.1109/focs.2017.97.
  * [8] R.F. Werner.  Large deviations and mean-field quantum systems.  In Quantum Probability and Related Topics, pages 349–381. WORLD SCIENTIFIC, jul 1992.  doi:10.1142/9789814354783_0024.
  * [9] Christina V. Kraus, Maciej Lewenstein, and J. Ignacio Cirac.  Ground states of fermionic lattice hamiltonians with permutation symmetry.  Physical Review A, 88(2), aug 2013.  doi:10.1103/physreva.88.022335.
  * [10] Fernando G. S. L. Brandão and Aram W. Harrow.  Product-state approximations to quantum states.  Communications in Mathematical Physics, 342(1):47–80, jan 2016\.  doi:10.1007/s00220-016-2575-1.
  * [11] Scott Aaronson and Lijie Chen.  Complexity-theoretic foundations of quantum supremacy experiments.  In 32nd Computational Complexity Conference (CCC 2017). Schloss Dagstuhl-Leibniz-Zentrum fuer Informatik, 2017.  arXiv:1612.05903.
  * [12] Madan Lal Mehta.  Random matrices, volume 142.  Elsevier, 2004.
  * [13] James Demmel, Ioana Dumitriu, and Olga Holtz.  Fast linear algebra is stable.  Numerische Mathematik, 108(1):59–91, oct 2007.  doi:10.1007/s00211-007-0114-x.
  * [14] Guang Hao Low and Isaac L. Chuang.  Optimal Hamiltonian simulation by quantum signal processing.  Phys. Rev. Lett., 118:010501, 2017.  arXiv:1606.02685v2, doi:10.1103/PhysRevLett.118.010501.
  * [15] Guang Hao Low and Isaac L Chuang.  Hamiltonian simulation by qubitization.  2016\.  arXiv:1610.06546.
  * [16] Dominic W. Berry, Andrew M. Childs, Richard Cleve, Robin Kothari, and Rolando D. Somma.  Exponential improvement in precision for simulating sparse Hamiltonians.  pages 283–292, 2014.  arXiv:1312.1414, doi:10.1145/2591796.2591854.
  * [17] Dominic W. Berry, Andrew M. Childs, Richard Cleve, Robin Kothari, and Rolando D. Somma.  Simulating Hamiltonian dynamics with a truncated Taylor series.  Phys. Rev. Lett., 114:090502, 2015.  arXiv:1412.4687, doi:10.1103/PhysRevLett.114.090502.
  * [18] D. W. Berry, A. M. Childs, and R. Kothari.  Hamiltonian simulation with nearly optimal dependence on all parameters.  In 2015 IEEE 56th Annual Symposium on Foundations of Computer Science, pages 792–809, Oct 2015.  arXiv:1501.01715, doi:10.1109/FOCS.2015.54.
  * [19] Alexei Yu Kitaev, Alexander Shen, and Mikhail N Vyalyi.  Classical and quantum computation, volume 47.  American Mathematical Society Providence, 2002.  doi:10.1090/gsm/047.
  * [20] Gilles Brassard, Peter Hoyer, Michele Mosca, and Alain Tapp.  Quantum amplitude amplification and estimation.  Contemporary Mathematics, 305:53–74, 2002.  doi:10.1090/conm/305.
  * [21] Anthony Carbery and James Wright.  Distributional and $l^$\\{$q$\\}$$ norm inequalities for polynomials over convex bodies in $$\\{$\bbb r$\\}$^n$.  Mathematical Research Letters, 8(3):233–248, 2001.  doi:10.4310/mrl.2001.v8.n3.a1.
  * [22] Dave Wecker, Matthew B. Hastings, Nathan Wiebe, Bryan K. Clark, Chetan Nayak, and Matthias Troyer.  Solving strongly correlated electron models on a quantum computer.  Physical Review A, 92(6), dec 2015.  doi:10.1103/physreva.92.062318.
  * [23] Matthew B. Hastings.  The asymptotics of quantum max-flow min-cut.  Communications in Mathematical Physics, 351(1):387–418, nov 2016\.  doi:10.1007/s00220-016-2791-8.

##  Appendix A Contracting Tensor Networks

In this appendix, we briefly explain how to estimate the average of a tensor
network over tensors with entries chosen from a Gaussian distribution, either
real or complex, in both cases with zero mean and unit variance. This is a
fairly standard technique, and Ref. [23] is an example that uses this
technique heavily and can be used as a reference for some of the ideas.

A tensor network is defined by some graph, with each vertex being labelled by
some tensor. The degree of each vertex is equal to the order of the
corresponding tensor. If the tensors are not symmetric, then for each vertex
we also specify some ordering on the edges attached to that vertex.

A tensor network $L$ with $m$ tensors evaluates to a scalar ${\rm Val}(L)$
that we will call the value of the tensor network. The value is computed as
follows: for each way of assigning an index to each edge (there are $N^{pm/2}$
such choices if there are $m$ tensors in the network each of order $p$), for
every vertex $v$ we compute some scalar; this scalar is equal to a particular
entry of the tensor $T^{v}$ corresponding to that vertex $v$. The entry that
is chosen is given by the indices on the edges attached to that vertex (if the
tensor is not symmetric, one uses the ordering on edges of that vertex to
determine the order of indices in the tensor). Then, the product of this
scalar over all vertices $v$ gives some scalar for each way of assigning an
index to each edge. Finally, the sum of this scalar over all possible
assignments to edges is ${\rm Val}(L)$, so that

${\rm Val}(L)=\sum_{\rm assignments}\prod_{{\rm vertex}\,v}(T^{\rm v})_{\rm
edges},$ (62)

where the sum is over choices of indices, the product is over vertices $v$,
and $(T^{v})_{\rm edges}$ denotes the entry of a tensor corresponding to the
given vertex $v$ determined by the indices on edges attached to that vertex.

Thus given $m$ tensors of order $p$, the value is an $m$-th order polynomial
in the entries. Equivalently, it is a sum over $mp$ different indices (one for
each leg attached to each tensor), with $mp/2$ $\delta$-functions (one for
each edge) constraining any pair of legs on two different tensors which are
connected by an edge to have the same index.

Consider a tensor network in which vertices are labelled by $G$ or
$\overline{G}$ in the complex ensemble or by $G$ in the real ensemble. We wish
to average the value of this network over Gaussian choices of $G$. Let us
first consider the case in which the tensors are not symmetrized. By Wick’s
theorem, this is computed by a sum over pairings. Each pairing is either (if
all tensors are chosen from a real Gaussian distribution) a way of
partitioning the set of tensors into distinct pairs or (if tensors are chosen
from the complex distribution) a way of pairing the set into distinct pairs
such that we only pair $G$ with $\overline{G}$. If there are $m$ tensors in
the network, then in the real case there are $m!!$ pairings and in the complex
case there are $(m/2)!$ pairings assuming that there are $m/2$ tensors $G$ and
$m/2$ tensors $\overline{G}$ (if not, the expectation value vanishes).

For each pairing, we compute a scalar as follows. For every pair $T$ and
$T^{\prime}$, replace the two tensors $T_{\mu_{1},\ldots,\mu_{p}}$ and
$T^{\prime}_{\nu_{1},\ldots,\nu_{p}}$ in the polynomial with a product of
$\delta$-functions $\prod_{a=1}^{p}\delta_{\mu_{a},\nu_{a}}$. After making
this replacement, sum over all indices. The result will be equal to
$N^{n_{loop}}$ where $n_{loop}$ is some positive integer. Summing
$N^{n_{loop}}$ over all pairings computes the expectation value of the value
of the network.

This can conveniently be represented graphically as follows: given a tensor
network, and given any tensors that are paired, “cut” those two tensors out of
the network, leaving $2p$ open edges in the network (the edges that were
attached to those tensors). Then, attach those edges to each other, attaching
each of the $p$ edges that were attached to the first tensor to the
corresponding one that was attached to the second tensor. After doing this for
every pair, what is left is some collection of $n_{loop}$ closed loops.

If the tensors are symmetrized, for each pair we must sum over $p!$ different
ways of attaching the edges that were attached to the first tensor to those
that were attached to the second tensor.

One can easily also consider the case in which we insert tensors $v_{\rm
sig}^{\otimes p}$ into the network. This can be done by applying the rotation
to make $v_{\rm sig}=(\sqrt{N},0,0,\ldots)$. Then, replace each tensor $v_{\rm
sig}^{\otimes p}$ with $\sqrt{N}$ multiplying a product of $\delta$-functions
on its edge, $\prod_{a=1}^{p}\delta_{\mu_{a},1}$. We can represent this
graphically simply by writing those $\delta$-functions as open edges; again
one sums $N^{n_{loop}}$ over pairings to compute the expectation value of the
value of the network, but any open edge with open ends does not count as a
closed loop.

##  Appendix B Second-Quantized Formalism

In this section, we briefly review the second-quantized formalism for bosons.
This formalism uses a countably infinite basis of states. In this appendix, we
write these basis states as $|n_{1},n_{2},\ldots,n_{N}\rangle$, where each of
the $n_{\mu}$ for $\mu\in\\{1,\ldots,N\\}$ is a non-negative integer, called
an “occupation number”. We may think of this as labelling possible states of
the tensor product of $N$ different harmonic oscillators, i.e., each $n_{\mu}$
labels a basis state for the $\mu$-th harmonic oscillator, sometimes called
the $\mu$-th “mode”.

The creation operator $a^{\dagger}_{\mu}$ is defined by

$a^{\dagger}_{\mu}|n_{1},\ldots,n_{N}\rangle=\sqrt{n_{\mu}+1}|n_{1},\ldots,n_{\mu-1},n_{\mu}+1,n_{\mu+1},\ldots,n_{N}\rangle,$

i.e., it increases $n_{\mu}$ by one, leaving the other $n_{\nu}$ for
$\nu\neq\mu$ unchanged, and it multiplies by a scalar $\sqrt{n_{\mu}+1}$. The
annihilation operator $a_{\mu}$ is defined to be the Hermitian conjugate of
the creation operator and one may verify that

$a^{\dagger}_{\mu}a_{\mu}|n_{1},\ldots,n_{N}\rangle=n_{\mu}|n_{1},\ldots,n_{N}\rangle.$

The number operator $n\equiv\sum_{\mu}a^{\dagger}_{\mu}a_{\mu}$ obeys
$n|n_{1},\ldots,n_{N}\rangle=\sum_{\mu}n_{\mu}|n_{1},\ldots,n_{N}\rangle$. In
this paper we restrict to the subspace with $n={n_{bos}}$. This gives us a
finite basis of states; indeed, restricting to any finite value for $n$ will
give a finite basis.

Finally, the single-particle density matrix is defined by

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

The remarkable property of the second-quantized notation is how it describes
behavior in the symmetric subspace of the full Hilbert space. Any Hamiltonian
such as Eq. (5) which is symmetric under permutation (that is, permuting basis
vectors by replacing
$|\mu_{1}\rangle\otimes\ldots\otimes|\mu_{{n_{bos}}}\rangle$ with
$|\mu_{\pi(1)}\rangle\otimes\ldots\otimes|\mu_{\pi({n_{bos}})}\rangle$ where
$\pi$ is an arbitrary permutation on ${n_{bos}}$ elements), will preserve the
symmetric subspace, i.e., the subspace which is invariant under permutation.
An othonormal basis for this symmetric subspace can be given by symmetrizing
the basis vectors $|\mu_{1}\rangle\otimes\ldots\otimes|\mu_{{n_{bos}}}\rangle$
for $\mu_{1}\leq\mu_{2}\leq\ldots\leq\mu_{n}$; we impose this (arbitrary)
ordering on the $\mu$ to avoid redundancies. The symmetrization of this basis
vector is then equal to (using the notation of this appendix) a basis vector
of the form $|n_{1},n_{2},\ldots,n_{N}\rangle$ where for each
$\mu\in{1,\ldots,N}$, the number $n_{\mu}$ is equal to the number of
$a\in{1,\ldots,{n_{bos}}}$ such that $\mu_{a}=\mu$. This basis is also equal
to (up to normalization constant) the basis $a^{\dagger}_{\mu_{1}}\ldots
a^{\dagger}_{\mu_{{n_{bos}}}}|0\rangle$ used in the text. Then, the single
particle density matrix that we have defined above using creation and
annihilation operators is equal to (using the formalism of the full Hilbert
space) the reduced density matrix (i.e., the marginal) on the first tensor
factor in the full Hilbert space (indeed, by symmetry, it is also equal to the
marginal on any of the ${n_{bos}}$ tensor factors). Similarly, a symmetric
Hamiltonian such as Eq. (5) can be expressed compactly using creation and
annihilation operators as in Eq. (6).