Statistical Estimation in the Spiked Tensor Model via the Quantum Approximate Optimization Algorithm

Authors: leo  zhou
joao  basso
song  mei
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 20:21:35

Abstract: The quantum approximate optimization algorithm (QAOA) is a general-purpose algorithm for combinatorial optimization that has been a promising avenue for near-term quantum advantage. 
In this paper, we analyze the performance of the QAOA on the spiked tensor model, a statistical estimation problem that exhibits a large computational-statistical gap classically. 
We prove that the weak recovery threshold of $1$-step QAOA matches that of $1$-step tensor power iteration. Additional heuristic calculations suggest that the weak recovery threshold of $p$-step QAOA matches that of $p$-step tensor power iteration when $p$ is a fixed constant. This further implies that multi-step QAOA with tensor unfolding could achieve, but not surpass, the asymptotic classical computation threshold $\Theta(n^{(q-2)/4})$ for spiked $q$-tensors. 
Meanwhile, we characterize the asymptotic overlap distribution for $p$-step QAOA, discovering an intriguing sine-Gaussian law verified through simulations. For some $p$ and $q$, the QAOA has an effective recovery threshold that is a constant factor better than tensor power iteration.
Of independent interest, our proof techniques employ the Fourier transform to handle difficult combinatorial sums, a novel approach differing from prior QAOA analyses on spin-glass models without planted structure.

Main paper: 
 introduction: 
 We study statistical estimation in the spiked tensor model, where we observe a q-tensor Y ∈ R n q in n q dimensions given byY = (λ n /n q/2 ) • u ⊗q + (1/ √ n) • W ∈ R n q . (1.1)Here u ∼ Unif({+1, -1} n ) is some hidden signal, 1 and W ∈ (R n ) ⊗q is a noise tensor whose entries are i.i.d. standard Gaussian N (0, 1). The parameter λ n > 0 is the signal-to-noise ratio (SNR). The goal is to estimate u given only access to Y . That is, we seek an estimator û : (R n ) ⊗q → S n-1 ( √ n) achieving nontrivial overlap with the signal:lim inf n→∞ E[⟨û(Y ), u⟩ 2 /n 2 ] > 0. (1.2)This task is known as weak recovery in the spiked tensor model.The spiked tensor model is a famous problem because it exhibits a huge computational-statistical gap, referring to regimes of SNR where the statistical estimation problem is information-theoreotically solvable, but no efficient algorithm has been found. For example, it is known that the Bayes optimal estimator achieves non-trivial overlap with the signal u when λ n > λ IT for some constant threshold λ IT = Θ (1) , whereas the problem is information-theoretically impossible when λ n ≤ λ IT [1] . Furthermore, the maximum likelihood estimator also achieves non-trivial overlap with the signal when λ n > λ MLE for some λ MLE = Θ (1) . However, the best-known polynomial-time classical algorithms for computing a non-trivial estimator require a much higher SNR of λ n = Θ(n (q-2)/4 ). These include tensor power iteration, gradient descent, approximate message passing, and spectral methods with tensor unfolding [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] . Indeed, assuming the secret leakage planted clique conjecture, Ref. [13] proves an Ω(n (q-2)/4 ) lower bound on the SNR needed by any polynomial-time classical algorithm. See Fig. 1 for an illustration of the different SNR thresholds and Section 2.1 for more background.On the other hand, quantum algorithms are widely believed to have computational advantages over classical algorithms for many problem classes. In particular, we focus on the Quantum Approximate Optimization Algorithm (QAOA) [14] , a general-purpose quantum optimization algorithm that can be applied to optimize any objective function on bit-strings. The QAOA has received an enormous amount of attention in the quantum computing community for several reasons. First, the QAOA is simple and allows efficient implementation on near-term quantum hardware with many applications [15] [16] [17] [18] [19] . Additionally, the QAOA is computationally universal [20] , and its generalization can realize other powerful algorithms such as the quantum singular value transform [21] . Under common complexity-theoretic assumptions, no classical device can efficiently simulate the output distribution of the QAOA even at shallow depth [22, 23] . Furthermore, the QAOA is guaranteed to find optimal solutions when its number of steps (or depth) diverges [14] . Nevertheless, analyzing the asymptotic performance of QAOA remains challenging: classical simulation of the algorithm is limited to small problem dimension n, and analytical computations are often highly non-trivial [24] [25] [26] [27] [28] . Given the enormous Ω(n (q-2)/4 ) computational-statistical gap in the spiked tensor model (compared to e.g., the constant factor gap in spin-glass optimization [29] ), it is an interesting open question whether the QAOA, as a realistic quantum algorithm with asymptotic convergence guarantees, can provide any computational advantages.In this work, we investigate the performance of QAOA for the spiked tensor model. In particular, we choose the log-likelihood objective of spiked tensor C(z) = ⟨Y , z ⊗q ⟩/n (q-2)/2 . Its maximizer, the maximum likelihood estimator, achieves non-trivial overlap with the signal whenever λ n > λ MLE = Θ (1) . While the infinite-step QAOA could compute the maximizer, we are interested in the performance of QAOA when the depth is polynomial in the problem size, and hope that it can surpass the Θ(n (q-2)/4 ) classical threshold. Although some limitations of the QAOA are known for certain random optimization problems in the low-depth regime [26, [30] [31] [32] , these negative results do not apply to the spiked tensor model since they rely on either sparse connectivity or concentration, both of which are absent in the current setting. Here, as a first attempt to bridge the gap in understanding how well a popular quantum algorithm may perform on a classically hard statistical estimation problem, we study the asymptotic behavior of the QAOA on the spiked tensor model in the constant-depth regime, where we are able to obtain rigorous and analytical results.Our contribution. In this paper, we analyze the signal-to-noise ratio threshold of p-step QAOA for weak recovery in the spiked tensor model, in the regime of fixed p and n approaching infinity. For p = 1, we prove the weak recovery threshold is λ n = Θ(n (q-1)/2 ), matching that of 1-step tensor power iteration. For p > 1, heuristic calculations suggest the threshold is λ n = Θ(n (q-2+εp)/2 ) log λ n /log n 0 (q -2)/4 (q -2)/2 (q -1)/2Information-theoretic threshold Classical computational threshold. Threshold for spectral algorithm and multi-step QAOA with tensor unfolding Multi-step tensor power iteration and multi-step QAOA 1-step tensor power iteration and 1-step QAOA where ε p = (q -2)/[(q -1) p -1], again matching p-step tensor power iteration. Additionally, given an initialization vector with n c /n correlation to the signal for 1/2 < c < 1, we prove the weak recovery threshold for 1-step QAOA is λ n = Θ(n (1-c)(q-1) ), identical to 1-step tensor power iteration. These results indicate that constant-step QAOA has the same asymptotic recovery threshold as tensor power iteration in the spiked tensor model. Meanwhile, further heuristic analysis suggests that QAOA with tensor unfolding could achieve the classical computation threshold Θ(n (q-2)/4 ).Furthermore, we derive the asymptotic distribution of the overlap for p-step QAOA, revealing an intriguing sine-Gaussian law distinct from p-step tensor power iteration. Analyzing the second moment, we see that, for certain (p, q) pairs, the QAOA effectively has a recovery threshold that is a constant factor better than tensor power iteration. Since there are classical algorithms that achieve better recovery thresholds than power iteration, it remains an interesting open question whether quantum advantage over the state-of-the-art classical algorithms may be obtained at larger QAOA depths that grow with system size. To our current knowledge, our work is the first to obtain analytical results using the QAOA for a statistical inference problem.The proof of the sine-Gaussian distribution adopted novel techniques, including using discrete Fourier transforms and the central limit theorem to handle combinatorial summations. The Fourier transform technique also allows us to replace nonlinear polynomials in the exponents with dual variables, leaving linear exponents that become easy in combinatorial sums. These techniques are of independent interest and could be useful for analyzing the QAOA in other models.2 Background and related work 
spiked tensor model and prior algorithms: 
 The spiked tensor model (1.1) was first introduced as a statistical model for tensor principal component analysis in [2] , where it was studied with a spherical prior u ∈ S n-1 ( √ n). The information-theoretic threshold for weak recovery under this model with the spherical prior [3, 7, 8] and the Rademacher prior u ∈ {±1} n [1] are both λ n = Θ(1).Tensor power iteration. A well-studied classical algorithm for the spiked tensor model is tensor power iteration [2, 10, 33] . Starting from a uniform random initialization û0 ∼ Unif(S n-1 ), the k-th iteration is given by ûk , whereûk = √ nY [û ⊗(q-1) k-1 ]/ Y [û ⊗(q-1) k-1 ] 2 , k ≥ 1, û0 ∼ Unif(S n-1 ). (2.1)Here, Y [û ⊗(q-1) ] ∈ R n denotes contracting the order-q tensor Y ∈ R n q with the order-(q -1) tensor û⊗(q-1) ∈ R n q-1 . It is shown that with (log n) iterations, weak recovery is possible if the SNR satisfies λ n = Ω(n (q-2)/2 / polylog(n)) [10, 33] . However, tensor power iteration does not match the best-known classical algorithms. Furthermore, we remark that rounding the tensor power iteration to sign(û k ) ∈ {±1} n does not give a better threshold.Other classical algorithms and related results. [2] showed that the tensor power iteration and approximate message passing algorithms with random initialization can recover the signal provided λ n = Ω(n (q-1)/2 ). This SNR threshold was later improved to λ n = Ω(n (q-2)/2 ) by [3, 10, 33] for these same methods. The same threshold λ n = Ω(n (q-2)/2 ) could also be achieved by gradient descent and Langevin dynamics as proved in [9] . On maximum likelihood estimation for the spiked tensor model with a spherical prior, [5, 6] studied the loss landscape, providing intuition that it contains many saddle points and local minima near the equator, but no bad critical points off the equator.The best currently known polynomial-time algorithms can achieve a sharp threshold of λ n = Ω(n (q-2)/4 ). These include spectral methods with tensor unfolding [2, 11] , sum-of-squares algorithms [34] [35] [36] , sophisticated iteration algorithms [37] [38] [39] , and gradient descent on the smoothed landscape [40, 41] .Another line of research has attempted to prove computational lower bounds in restricted computational models, including low-degree polynomials and statistical query algorithms [42, 43] . Under the secreted leakage planted clique conjecture, [13] proved that any classical polynomial-time algorithm requires λ n = Ω(n (q-2)/4 ) for weak recovery of the signal.A quantum algorithm by Ref. [44] . To the best of our knowledge, the only prior quantum algorithm proposed for the spiked tensor model with provable guarantees is by Hastings in Ref. [44] . Hastings' algorithm is based on a spectral method for a Hamiltonian on M bosons over n modes, living in a Hilbert space of dimension n M , where M ≫ [n (q-2)/4 /λ n ] 4/(q-2) × polylog(n). Finding the dominant eigenvector of this Hamiltonian allows for weak recovery in the regime where λ n = Θ(n (q-2)/4 ). In this regime, where M = Ω(polylog n), the standard classical matrix power iteration algorithm can extract the dominant eigenvector and recover the signal in Õ(n M ) time. For the proposed quantum algorithm, Ref. [44] uses a combination of quantum phase estimation, amplitude amplification, and clever state initialization to recover the signal in Õ(n M/4 ) time, achieving a quartic speedup. (A few months after our paper appeared online, a related work [45] emerged, simplifying Hastings' algorithm and generalizing it to another planted inference problem.)We remark that Hastings' algorithm runs in superpolynomial time n Ω(polylog n) and does not improve over the asymptotic computational threshold in SNR for recovery (although a constant factor improvement is possible). For comparison, the classical spectral method based on tensor unfolding [2, 11] achieves recovery when λ n > n (q-2)/4 in polynomial time O(poly(n q )). In this work, we study the QAOA in the constant-step regime, where the gate complexity grows only linearly in the problem size O(n q ).The quantum approximate optimization algorithm (QAOA) was introduced by [14] as a quantum algorithm for finding approximate solutions to combinatorial optimization problems. The QAOA can be applied to optimize any cost function on bit-strings, C : {±1} n → R. In the spiked tensor model, we consider optimizing the log-likelihood function given byûMLE = arg max σ∈{±1} n C(σ) = ⟨Y , σ ⊗q ⟩/n (q-2)/2 . (2.2)The maximum likelihood estimator ûMLE achieves non-trivial correlation with the signal when λ n > λ MLE for some constant λ MLE = Θ(1). However, classical algorithms cannot efficiently compute the MLE unless λ n = Ω(n (q-2)/4 ) [13] . This paper investigates whether QAOA could compute ûMLE , or an approximate estimator, for smaller values of λ n .The inputs to the QAOA algorithm are a cost function C :{±1} n → R and parameter vectors γ, β ∈ R p . The initial QAOA state |s⟩ = 2 -n/2 z |z⟩ is the rescaled all-one vector 2 -n/2 1 2 n ∈ C 2 n, assigning equal probability to measuring each possible bit-string upon quantum measurement. See Appendix A.1 for a review of quantum computing terminology, where we also define the Pauli operators {X k , Y k , Z k } n k=1 acting on the k-th qubit. The cost function C associates with a 2 n × 2 n diagonal matrix, where the |z|'th diagonal gives C(z). For the spiked tensor model with cost functionC(z) = ⟨Y , z ⊗q ⟩/n (q-2)/2 , this matrix is C = n j1,...,jq=1 Y j1•••jq Z 1 • • • Z q /n (q-2)/2 ∈ C 2 n ×2 n . Letting B = n j=1 X j ∈ C 2 n ×2 n, for any parameter (γ, β), the unitary matrices e -iγC , e -iγB ∈ C 2 n ×2 n are matrix exponents of -iγC and -iγB. Given γ, β ∈ R p , the p-step QAOA state is|γ, β⟩ = e -iβpB e -iγpC • • • e -iβ1B e -iγ1C |s⟩ ∈ C 2 n . (2.3) One can verify |γ, β⟩ is a unit vector since |s⟩ ∈ C 2 n is unit and e -iβ k B ∈ C 2 n ×2 n and e -iγ k C ∈ C 2 n ×2 nare unitary matrices. After preparing the quantum state |γ, β⟩, QAOA samples a bit string z ∼ |γ, β⟩ in {±1} n by quantum measurement. In our main results, we will analyze the distribution of the overlap R QAOA of this quantum measurement z with respect to the signal u: R QAOA ≡ z ⊤ u/n = 1 n n i=1 z i u i ∈ [-1, 1]. (2.4) For any function f (z) = n k=0 (j1,••• ,j k ) fj1•••j k z j1 • • • z j k ,(Z) = n k=0 (j1,••• ,j k ) fj1•••j k Z j1 • • • Z j k ∈ R 2 n ×2 nfor Pauli-Z matrices Z j . To simplify the notations, we denote ⟨•⟩ γ,β by the expectation with the quantum measurement from |γ, β⟩, so that ⟨f (Z)⟩ γ,β = ⟨γ, β|f (Z)|γ, β⟩.(2.5) In the main theorems of this paper, we will focus on the second moment of the overlap of QAOA, denoted as⟨R 2 QAOA ⟩ γ,β = ⟨γ, β| R 2 |γ, β⟩, where R ≡ 1 n n i=1 u i Z i .We defer further related literature on theoretical analyses of the QAOA to Appendix A.2.In terms of experimental realizations, the QAOA has been implemented in quantum computing platforms such as trapped ions [16, 19] , superconducting qubits [17] , and neutral atoms [18] , for optimization problems with up to 179 bit variables. Implementing the QAOA for the spiked tensor model, however, poses additional challenges due to the all-to-all connectivity in its cost function, leading to a higher overhead in the number of quantum gates and circuit compilation costs. Currently, the largest experimental implementations for problems with dense connectivity include 17-bit Sherrington-Kirkpatrick spin-glass models on superconducting qubits [17] , and an 18-bit LABS problem on trapped-ion quantum processors [19] . We expect larger problems can be implemented as quantum hardware matures, but quantum error-correction is likely necessary to observe any quantum advantage at scale [46] . (γ n , β n ) ∈ R >0 × [0, 2π].The quantum state |γ n , β n ⟩ depends randomly on Y through C(σ) = ⟨Y , σ ⊗q ⟩/n (q-2)/2 . Our main results characterize the distribution of the overlap R QAOA = û⊤ u/n between a sample û ∼ |γ n , β n ⟩ and the signal vector u. Theorem 1 (Weak recovery threshold and overlap distribution for 1-step QAOA). Consider the spiked tensor model (1.1) and the 1-step QAOA overlap as defined above. Then the following hold.(a) Take any sequence of {γ n } n≥1 ⊆ R, {β n } n≥1 ⊆ [0, 2π], and any sequence of {λ n } n≥1 ⊆ [0, ∞) with lim n→∞ λ n /n (q-1)/2 = 0. We haveEQUATION(b) Take any sequence of {γ n } n≥1 , {β n } n≥1 , and {λ n } n≥1 which satisfieslim n→∞ (γ n , β n , λ n /n (q-1)/2 ) = (γ, β, Λ). (3.2)Then, over the randomness of Y and the quantum measurement, the overlap R QAOA of the 1-step QAOA converges in distribution to a sine-Gaussian law asR QAOA d-→ e -2qγ 2 sin(2β) sin(2qΛγG q-1 ), where G ∼ N (0, 1).(c) As a corollary of (b), under the asymptotic limit of (3.2) with Λ > 0, γ > 0, and β ̸ ∈ {kπ/2 : k ∈ Z}, we havelim n→∞ E Y [⟨R 2 QAOA ⟩ γn,βn ] > 0. (3.4)The full proof of Theorem 1 is contained in Appendix C.Remark 3.1 (Weak recovery threshold). Theorem 1(c) implies that when λ n = Θ(n (q-1)/2 ), the overlap will be non-zero with non-trivial probability over both the random draw of the tensor and the quantum randomness. In contrast, Theorem 1(a) shows that when λ n = o(n (q-1)/2 ) the overlap will be zero with high probability. This establishes that λ n = Θ(n (q-1)/2 ) is the weak recovery threshold of 1-step QAOA in the spiked tensor model. Remark 3.2 (Overlap distribution). Theorem 1 does not show that the overlap distribution for a typical instance Y converges to the same sine-Gaussian law. In Section 4, we perform numerical simulations that provide evidence that the overlap distribution will concentrate over the random draw of Y , which would imply that the overlap distribution is indeed sine-Gaussian for any typical Y .Comparison with classical tensor power iteration. The 1-step tensor power iteration estimator (Eq. (2.1)) is redefined here for the reader's convenience:û1 = √ nY [û ⊗(q-1) 0 ]/∥Y [û ⊗(q-1) 0 ]∥ 2 , where û0 ∼ Unif(S n-1) is a random initialization vector. In the following proposition, we show that the weak recovery threshold for the 1-step power iteration estimator is also λ n = Θ(n (q-1)/2 ), and we provide the distribution of the overlap R PI ≡ û⊤ 1 u/n between the power iteration estimator û1 and the signal u. Proposition 3.3 (Weak recovery threshold for 1-step tensor power iteration). Assume that the rescaled signal-to-noise ratio has a limit lim n→∞ λ n /n (q-1)/2 = Λ. Then over the randomness of W and initialization û0 , the overlap R PI of the power iteration estimator with the signal converges in distribution to R PI d -→ sin[arctan(ΛG q-1 )], where G ∼ N (0, 1).(3.5)As a corollary, when lim n→∞ λ n /n (q-1)/2 = 0, we have R PI p -→ 0.The proof of Proposition 3.3 is contained in Appendix H.1. Remark 3.4 (Comparing the overlaps). Theorem 1 and Proposition 3.3 show that both 1-step QAOA and 1-step power iteration have the same weak recovery threshold λ n = Θ(n (q-1)/2 ). To compare the two algorithms more precisely, we take the limit lim n→∞ λ n /n (q-1)/2 = Λ for some small Λ > 0. Eq. (3.3) and Eq. (3.5) give the limiting squared overlap distributions for 1-step QAOA and 1-step power iteration, respectively:EQUATIONThis givesEQUATIONwhere the maximizer is(γ ⋆ , β ⋆ ) = ( 1 2√ q , π/4). Thus, for q > e, 1-step QAOA gives better overlap than 1-step power iteration. Remark 3.5 (Rounding via sign(û) will not improve the overlap). The readers may wonder whether the overlap of tensor power iteration will be improved by rounding the estimator via ū1 = sign(û 1 ) ∈ {±1} n , outputting an estimator in the signal space. Defining R PI = ū⊤ 1 u/n, it is straightforward to show that as lim n→∞ λ n /n (q-1)/2 = Λ,R PI d -→ Φ(ΛG q-1 ), where G ∼ N (0, 1), Φ(t) = 2 × P Z∼N (0,1) (Z ≤ t) -1. (3.8)Hence, the computational threshold has the same exponent by rounding, but the overlap becomes smaller:EQUATIONRemark 3.6 (Sine-Gaussian law versus sine-arctan-Gaussian law). The sine-Gaussian law of QAOA is particularly interesting in that the overlap will not concentrate as Λ → ∞. Instead, it will satisfy a sine-uniform distribution, i.e., sin(2qΛγG q-1 ) d -→ sin(U ) for U ∼ Unif([0, 2π]). In contrast, the sine-arctan-Gaussian law of tensor power iteration will concentrate at {±1} as Λ → ∞.In Appendix E, we also study the scenario where prior information about the signal may be leveraged to recover the signal with a smaller SNR. There, we rigorously analyzed the 1-step QAOA applied to boost the signal in a weak estimator in Theorem 2, and compared it classical power iteration in Proposition E.2. Our result shows that the 1-step QAOA has the same asymptotic computational efficiency as 1-step power iteration, albeit with a constant-factor better overlap in the Λ ≪ 1 regime when q > e.3.2 Weak recovery threshold and overlap distribution for p-step QAOA We next consider the general p-step QAOA for weak recovery in the spiked tensor model. Although it is known that the QAOA is able to output the MLE that weakly recovers the signal when p grows unboundedly with n, here we focus on a more analytically tractable regime where p is an arbitrary fixed constant in the n → ∞ limit. Using a physics-style derivation, we show that the p-step QAOA can achieve weak recovery when the signal-to-noise ratio satisfies EQUATIONEQUATIONThen, there are parameter-dependent coefficients (a p (γ, β), b p (γ, β)) such that over the randomness of Y and the quantum measurement, the overlap R QAOA of the p-step QAOA converges in distribution to a sine-Gaussian law asR QAOA d -→ a p sin(b p Λ 1/εp G (q-1) p ),where G ∼ N (0, 1).(3.12)The derivation of Claim 3.7 is contained in Appendix D. We remark that our derivation uses nonrigorous heuristics from physics such as the Dirac delta function and its Fourier transform to linearize exponents in combinatorial sums (see Appendix D.1 for a sketch). Analytical expressions for the coefficients a p (γ, β) and b p (γ, β) can be found in Appendix D.5.Remark 3.8 (Weak recovery threshold). As Λ → 0, Eq. (3.12) implies that R QAOA p -→ 0. Thus, Claim 3.7 implies that λ n = Θ(n (q-2+εp)/2 ) is the weak recovery threshold by the p-step QAOA in the spiked tensor model in the regime of fixed QAOA parameter. We believe this scaling is also the weak recovery threshold for the QAOA with any sequence of parameters (γ n , β n ), but proving this requires ruling out better performance of the QAOA when (γ n , β n ) is allowed to depend strongly on n as we have done in Theorem 1(a); we leave this as future work. Since ε p → 0 as p → ∞, this means λ n = Θ(n (q-2)/2 ) is the recovery threshold given a diverging number of QAOA steps (but constant with respect to n). However, this does not achieve the Θ(n (q-2)/4 ) computational threshold for classical algorithms.Comparison with classical tensor power iteration. We now compare the overlap from the p-step QAOA to that from the classical p-step tensor power iteration algorithm. We show that the weak recovery threshold for the p-step power iteration estimator is also λ n = Θ(n (q-2+εp)/2 ), and we provide the distribution of the overlap R PI ≡ û⊤ p u/n between the p-step power iteration estimator ûp (see Eq. (2.1)) and the signal u. Proposition 3.9 (Corollary of Lemma 3.2 of [33] ). Consider a random instance of the spiked tensor model with lim n→∞ λ n /n (q-2+εp)/2 = Λ. The overlap R PI of the p-step tensor power iteration algorithm converges in distribution asR PI d -→ sin arctan(Λ 1/εp G (q-1) p) , where G ∼ N (0, 1).(3.13)The proof of Proposition 3.9 is contained in Appendix H.3. Remark 3.10 (Comparing the overlaps). In the small Λ ≪ 1 regime, we haveR QAOA ≍ (|a p b p | εp Λ) 1/εp G (q-1) p and R PI ≍ Λ 1/εp G (q-1) p . (3.14)When |a p b p | > 1, the QAOA has a constant factor advantage over the classical power iteration algorithm in the overlap achieved, assuming the conjectured Claim 3.7 based on heuristic derivations.To quantify this advantage, we consider the quantum enhancement factor, |a p b p | εp , which is the factor that the signal-to-noise ratio can shrink for the QAOA while maintaining the same overlap as the power iteration algorithm. Effectively, this factor |a p b p | εp corresponds to a quantum improvement in the recovery threshold by the QAOA over classical power iteration. We numerically optimize |a p b p | εp with respect to the QAOA parameters (γ, β), and present the optimized values in Table 1 .Table 1 : The quantum enhancement factor |a p b p | εp of the p-step QAOA over the p-step tensor power iteration, for spiked q-tensors when λ n = Λn (q-2+εp)/2 in the Λ ≪ 1 regime. Note in the first row, which corresponds to p = 1 with ε 1 = 1, we know the optimal value |a 1 b 1 | = q/e from Eq. (3.7). The remaining values are optimized via a quasi-Newton method starting with 1000 heuristic initial guesses of (γ, β) and keeping the best value; hence, they currently should be considered as lower bounds on the best possible enhancement factors. Although neither the constant-step QAOA nor the tensor power iteration matches the Θ(n (q-2)/4 ) recovery threshold for the best polynomial-time classical algorithms, we can achieve this threshold using the idea of tensor unfolding. When q is even, the tensor Y ∈ R n q can be unfolded into a matrix Y :pY = (λ n /n q/2 ) • ūū ⊤ + (1/ √ n) • W ∈ R n q/2 ×n q/2 . (3.15)Here Y (j1,...,j q/2 ),(j q/2+1 ,...,jq) = Y j1•••jq , W (j1,...,j q/2 ),(j q/2+1 ,...,jq) = W j1•••jq , and ū = vec(u ⊗(q/2) ) ∈ {±1} n q/2 . Existing work [2, 47] have demonstrated that the leading eigenvector z of Y has non-vanishing correlation with the signal ū as soon as λ n > n (q-2)/4 . Furthermore, for the eigenvector z in such a regime, standard analysis as in [2] implies that the top singular vector of mat(z) ∈ R n×n q/2-1will have non-trivial overlap with the signal u, achieving the Θ(n (q-2)/4 ) weak recovery threshold for the spectral method with tensor-unfolding.A similar tensor-unfolding pre-processing could be applied to the QAOA to improve the computational threshold. Indeed, the QAOA method could be adopted to maximize the cost functionC( σ) = σ⊤ Y σ/n (q-1)/2 with decision variable σ ∈ {±1} n q/2. Notice that such a QAOA method needs to be applied to a n q/2 -qubit system. Effectively, C( σ) could be interpreted as the cost function of a spiked 2-tensor model of size n = n q/2 and with a rescaled signal-to-noise ratio λn = λ n /n (q-2)/4 . According to Claim 3.7, p-step QAOA outputs a long bit-string z ∈ {±1} n q/2 overlapping with the signal ū as long as λn = nεp/2 for ε p = 1/p. Translating to the scaling of λ n , the computational threshold for QAOA with tensor unfolding is λ n = Ω(n (q-2+ε ′ p )/4 ) where ε ′ p = q/p. This recovers the classical Θ(n (q-2)/4 ) threshold as p → ∞. 
numerical simulations: 
 We now validate our theoretical results by conducting numerical simulations of the QAOA through classical computers. In this section, we focus on the case of 1-step QAOA (p = 1) for the spiked matrix model (q = 2), where we can obtain an explicit formula the expected squared overlap at any finite problem dimension n (see Appendix F for a derivation):E Y [⟨R 2 QAOA ⟩ γ,β ] = n -1 2n e -8γ 2 (n-2)/n sin 2 (2β)[1 -cos n-2 (8λγ/n)] + n -1 n e -4γ 2 (n-1)/n sin(4β) sin(4λγ/n) cos n-2 (4λγ/n) + 1 n . (4.1)In Fig. 2 (a), we report the overlap distribution of 1-step QAOA (p = 1) for the spiked matrix model (q = 2) where the SNR is chosen as The top row shows data from 40 random 26-bit instances with q = 2 and λ n = n 1/(2p) . The bottom row shows data from 40 random 23-bit instances with q = 3 and λ n = n [1+1/(2 p -1)]/2 . Different columns correspond to different p, using the QAOA parameters (γ, β) that optimized |a p b p | εp in Table 1 . Dash gray lines connect data from the same instance. Blue histograms are the theoretical sine-Gaussian distributions in the n → ∞ limit, where R QAOA ∼ a p sin[b p G (q-1) p ] according to Claim 3.7. (Note here Λ = 1.) that, despite some finite sample effects, the predicted sine-Gaussian distribution matches the QAOA simulation.λ n = n 1/2 .Fig. 2 (b) reports the expected squared overlap from the QAOA simulations. The green dashed line is the theoretical prediction in the n → ∞ limit. The blue solid line is the finite n theoretical prediction from Eq. (4.1). The gray dots are the squared overlaps from individual QAOA instances simulated classically. The average over instances (red crosses) agrees well with the finite n theory prediction, which converges to the n → ∞ limit with order 1/n deviation.We also perform simulations for 1 ≤ p ≤ 5 and q = 2, 3. Fig. 3 plots the overlap distribution for p-step QAOA. The simulation curves follow the shape of the theoretical histograms for p ≤ 2. For p ≥ 3, the shapes of the simulated and theoretical overlap distributions do not match well, likely due to finite size effects (simulations for large n > 26 are computationally challenging).In Appendix G, we present additional numerical simulation results on higher (p, q) and find the second moment of the QAOA overlap converges to our theoretical predictions up to O(1/n) deviations. We also describe more details of the simulation methods.An interesting phenomenon apparent from Fig. 2 (a) and Fig. 3 is that the output distribution of the QAOA appears to concentrate over the randomness of instances Y , but not over the quantum measurements. This is in stark contrast to previous concentration results on the QAOA where concentration over measurements were shown, e.g., for spin-glass models in [24, 26, 32] . We note that such anti-concentration is also expected in the limit of zero noise (λ → ∞), where it is known the constant-p QAOA can prepare the GHZ state [48] . Since existing limitations of both classical [29] and quantum algorithms [26, [30] [31] [32] on various problems over random structures rely heavily on concentration, extending these negative results to the QAOA for the spiked tensor model do not seem possible due to the absence of concentration. Nevertheless, our analysis shows that the constant-p QAOA is unable to improve the recovery threshold in the spiked tensor model achieved by classical algorithms by more than a constant factor. 
discussion: 
 In this paper, we have investigated the power of quantum algorithms for the spiked tensor model, a canonical problem in statistical inference with a large computational-statistical gap that has so far eluded classical algorithms. We gave the first rigorous study of a polynomial-time quantum algorithm on this problem by analyzing the performance of the QAOA, a popular variational quantum algorithm that has been implemented on current quantum computing hardware. We showed that p-step QAOA achieves the same asymptotic SNR threshold for weak recovery as p-step tensor power iteration. A heuristic analysis showed that multi-step QAOA with tensor unfolding could achieve, but not surpass, the classical computation threshold Θ(n (q-2)/4 ). This implies that achieving a strong quantum advantage via the QAOA requires using a number of steps p that grows with n. However, we revealed that the asymptotic overlap distribution of QAOA exhibits an intriguing sine-Gaussian law, distinct from tensor power iteration. For certain parameters (p, q), the QAOA effectively has a recovery threshold that is a constant factor better, indicating a modest quantum advantage over the classical power iteration. Overall, while achieving identical scalings as power iteration, the QAOA demonstrates qualitative differences and potential for quantum speedups.There are many interesting questions that remain open. One worthy challenge would be a rigorous proof for the p > 1 analysis without relying on heuristic arguments. Additionally, it would be interesting to prove that the sine-Gaussian distribution is concentrated over problem instances but not over measurements, as suggested by our simulations. This is in contrast to recent results showing that the low-depth QAOA is concentrated over measurements [26, 32] , a seemingly essential ingredient for many proofs of algorithmic limitations [26, [29] [30] [31] [32] . Despite the absence of concentration in the spiked tensor setting, our results show that the constant-p QAOA has limited power, similar to the message of recent works [26, [30] [31] [32] [49] [50] [51] proving limitations up to p = O(log n). This suggests that demonstrating strong quantum advantage requires analyzing super-logarithmic depth QAOA, which remains an outstanding open question. Finally, it would be interesting to study quantum algorithms in other statistical inference models that classically exhibit computational-statistical gaps, including planted clique, Bayesian linear models, and sparse PCA. Overcoming any such gap with a polynomial-time quantum algorithm would be an exciting superpolynomial quantum speedup with practical relevance. n -dimensional unit complex vector ψ ∈ C 2 n satisfying i∈[2 n ] |ψ i | 2 = 1. Each bit-string z ∈ {±1} n associates with a quantum state |z⟩ ∈ C 2 n , representing the |z|'th canonical basis vector [0, • • • , 0, 1, 0, • • • , 0] ⊤ ∈ C 2 n, where only position |z| equals 1 (with |z| = 1 + j∈[n] 2 j-1 (1z j )/2 denoting the rank of bit-string z). Therefore, ψ = z∈{±1} ψ |z| |z⟩ where |ψ |z| | 2 gives the probability of observing z upon measurement. This represents ψ as a probability distribution over all 2 n bit-strings in {±1} n . The Pauli operators σ x , σ y , σ z on a single qubit are represented as 2 × 2 complex matrices:I = 1 0 0 1 , σ x = 0 1 1 0 , σ y = 0 -i i 0 , σ z = 1 0 0 -1 . (A.1)In an n-qubit system, the Pauli operators{X k , Y k , Z k } ∈ C 2 n ×2 n associated to the k-th qubit are defined by I ⊗(k-1) ⊗ {σ x , σ y , σ z } ⊗ I ⊗(n-k) ∈ C 2 n ×2 n, where ⊗ is the Kronecker product operator.A.2 Related literature on the QAOAIn terms of theoretical analysis of its computational complexity, the performance of the QAOA has been studied for various models, including the Sherrington-Kirkpatrick model [24] , MaxCut [25, 27] , the Max-q-XORSAT for regular hypergraphs [25] , q-spin spin-glass models [26, 52] , the ferromagnetic Ising model [53] , and random constraint satisfaction problems [28] . While [24] shows promising evidence for the QAOA to achieve the ground state energy of the Sherrington-Kirkpatrick model [54] , [26] proves that constant-step QAOA cannot achieve the ground state for q-spin spin-glass models in general.There is also a line of work aiming to prove computational hardness results for the QAOA and related quantum algorithms. [30, 31, 50, 51] studied the limitation of local quantum algorithms like the QAOA for solving combinatorial optimization problems on sparse random graphs, using the bounded light-cone of the algorithms at sufficiently low depths. This limitation was translated to the dense spin-glass models in [26] . Furthermore, [32, 49] proved hardness results for the QAOA by exploiting the symmetry of the problem. Note all previously known limitations of the QAOA in the average-case setting [26, [30] [31] [32] have relied on concentration of the output distribution in the Hamming weight basis, which is not present in the spiked tensor model.Our work studies the QAOA for a statistical inference problem, distinct from these existing results. Furthermore, we develop new techniques for analyzing the QAOA that do not exist in prior work.B Moment generating function of the QAOA overlap at p = 1We dedicate this section to derive a combinatorial expression for expected moment-generating function of the QAOA overlap, defined asM n (ζ; γ n , β n , λ n ) := ⟨e ζ R ⟩ γn,βn . (B.1)We write M n (ζ) = M n (ζ; γ n , β n , λ n ) for short. This quantity will be used in future derivations.We use the techniques and conventions first introduced in [24] . First, we define bistrings a ∈ B := {±1} 3 indexed as a = (a 1 , a m , a 2 ). Since p = 1, we write β = β 1 , γ = γ 1 . Additionally, define the quantities given byQ a = 1 2 ⟨a 1 |e iβX |1⟩ ⟨1|e -iβX |a 2 ⟩ , (B.2) Φ a = γ(a 1 -a 2 ). (B.3)We may also write {n a } a∈B ⊆ Z |B| where a∈A n a = n to assign a count to each bit-string. If we underscore the bit-string, we mean a = (a 1 , a 2 , . . . , a q ) ∈ B q . We also write Φ a = Φ a1a2•••aq .We now have the notation to state the following lemma. Lemma B.1 (QAOA overlap expected moment-generating function in the configuration basis for p = 1). The expectation over the spiked tensor disorder in Eq. (1.1) of the moment-generating function defined in Eq. (B.1) for p = 1 is given by E Y [M n (ζ)] = {na} n {n a } a∈B Q na a exp - 1 2n q-1 a∈B q Φ 2 a q s=1 n as + iλ n n q-1 a∈B q Φ a q s=1 (a s ) m n as + ζ n v∈B v m n v . (B.M n (ζ) = z 1 ,z m ,z 2 ⟨s| e iγC |z 1 ⟩ ⟨z 1 | e iβB e ζ R |z m ⟩ ⟨z m | e -iβB |z 2 ⟩ ⟨z 2 | e iγC |s⟩ = 1 2 n z 1 ,z m ,z 2 ⟨z 1 | e iβB |z m ⟩ e iγC(z 1 ) e ζ R(z m ) e iγC(z 2 ) ⟨z m | e -iβB |z 2 ⟩ = 1 2 n z 1 ,z m ,z 2 f * β (z 1 z m )f β (z m z 2 ) exp iγ(C(z 1 ) -C(z 2 )) + ζR(z m ) = 1 2 n z 1 ,z m ,z 2 f * β (z 1 z m )f β (z m z 2 ) × exp iγ n i1,...,iq=1 λ n n q-1 + W i1,...,iq n (q-1)/2 (z 1 i1 • • • z 1 iq -z 2 i1 • • • z 2 iq ) + ζ n n j=1 z m j , (B.6)where we defined f β (zz ′ ) = ⟨z| e -iβB |z ′ ⟩ since this quantity only depends on the bitwise product zz ′ . We also used the definitions of C(z) = ⟨Y , z ⊗q ⟩/n (q-2)/2 and R(z) = 1 n n j=1 z j . Next, we tranform the z j as follows:z 1 → z 1 z m , z 2 → z 2 z m . (B.7)This givesM n (ζ) = 1 2 n z 1 ,z m ,z 2 f * β (z 1 )f β (z 2 ) × exp iγ n i1,...,iq=1 λ n n q-1 + W i1,...,iq n (q-1)/2 z m i1 • • • z m iq (z 1 i1 • • • z 1 iq -z 2 i1 • • • z 2 iq ) + ζ n n j=1 z m j = 1 2 n z 1 ,z m ,z 2 f * β (z 1 )f β (z 2 ) × exp i n i1,...,iq=1 λ n n q-1 + W i1,...,iq n (q-1)/2 z m i1 • • • z m iq Φ i1,...,iq (Z) + ζ n n j=1 z m j , (B.8)where we denoted Z = (z 1 , z 2 ) andΦ i1,...,iq (Z) = γ(z 1 i1 • • • z 1 iq -z 2 i1 • • • z 2 iq ). (B.9)Hence, the expected moment-generating function isE Y [M n (ζ)] = 1 2 n z 1 ,z m ,z 2 f * β (z 1 )f β (z 2 ) × exp q i1,...,iq=1 iλ n n q-1 z m i1 • • • z m iq Φ i1,...,iq (Z) - 1 2n q-1 Φ 2 i1,...,iq (Z) + ζ n n j=1 z m j .(B.10)Now we change to the so-called configuration basis. For any bit-string 1 ≤ j ≤ n, we look at a new bit-string:(z 1 j , z m j , z 2 j ) ∈ B. (B.11)For any a ∈ B, we represent by n a the number of times that configuration a happens. Note that a∈B = n. For more details, we again refer the reader to [26, Appendix D.2] . Now, instead of counting over each bit of z 1 , z m , z 2 , we can count over configurations in A:E Y [M n (ζ)] = {na} n {n a } a∈B Q na a exp iλ n n q-1 a1,...,aq∈B Φ a1•••aq (a 1 ) m • • • (a q ) m n a1 • • • n aq - 1 2n q-1 a1,...,aq∈B Φ 2 a1•••aq n a1 • • • n aq + ζ n a∈B a m n a , (B.12)which finishes the proof of Lemma B.1.C.1 Proof sketch for Theorem 1(b) and emergence of sine-Gaussian law.Here we briefly sketch the proof of Theorem 1 for 1-step QAOA, explain how the sine-Gaussian law appears, and highlight the technical ideas.To derive the distribution of the QAOA overlap, we compute its expected moment-generating function.We start by following the steps from [24, 26] to reformulate the expected moment-generating function (MGF). With some algebra, we arrive at the following equation (see Appendix B and Lemma C.1 for the derivation):E Y [⟨e ζ R ⟩ γ,β ] = n t=0 n t sin(2β) 4n t e -γ 2 [n q -(n-2t) q ]/n q-1 S n,t (C.1) S n,t = n t ni=t t {n i } (-i) n1-n2+n3-n4 e (ζ/n)(n1+n2-n3-n4) Z n,t (n 1 -n 2 -n 3 + n 4 ), (C.2) Z n,t (k) = 1 2 n-t τ++τ-=n-t n -t τ + , τ - (e ζ/n cos 2 β + e -ζ/n sin 2 β) τ+ (e -ζ/n cos 2 β + e ζ/n sin 2 β) τ- × e iΛγ[((τ+-τ-)+k) q -((τ+-τ-)-k) q ]/n (q-1)/2 . (C.3)Looking upon the term Z n,t (k), we can interpret the summand τ + as a binomial variable Binom(nt, 1/2), and by the Central Limit Theorem, we have (τ+ -τ -)/ √ n d -→ G ∼ N (0, 1). This giveslim n→∞ Z n,t (k) = E G∼N (0,1) [e i2qkΛγG q-1 ] =: Z t (k).This is the step where the power of Gaussian appears. Next, assuming that we can replace Z n,t by Z t in the expression of S n,t as in (C.2), and using the multinomial theorem, we getS n,t • = E G∼N (0,1) n t ni=t t {n i } (-i) n1-n2+n3-n4 e (ζ/n)(n1+n2-n3-n4) e (n1-n2-n3+n4)i2qΛγG q-1 = E G∼N (0,1) {[4n sinh(ζ/n) sin(2qΛγG q-1 )] t } → E G∼N (0,1) {[4ζ sin(2qΛγG q-1 )] t } =: S t .This is the step where the sine-Gaussian distribution appears. Finally, suppose that we can replace S n,t by S t in (C.1), and using the Taylor expansion of the exponential function, we getE Y [⟨e ζ R ⟩ γ,β ] • = n t=0 n t sin(2β) 4n t e -γ 2 [n q -(n-2t) q ]/n q-1 E G∼N (0,1) [4ζ sin(2qΛγG q-1 )] t • → E ∞ t=0 1 t! [ζe -2qγ 2sin(2β) sin(2qΛγG q-1 )] t = E G∼N (0,1) e ζe -2qγ 2 sin(2β) sin(2qΛγG q-1 ) .This gives the moment-generating function of the sine-Gaussian law.We should notice that several steps in the above proof sketch are non-rigorous, in the sense that we could not sequentially take n → ∞ in Z n,t , S n,t , and the MGF. To make this step rigorous, we use the idea of discrete Fourier transform in Eq. (C.2) to decouple the two terms e (ζ/n)(n1+n2-n3-n4) and Z n,t (see Lemma C.1), which allows one to treat the n → ∞ limit of these two terms separately in the expression of S n,t . For more details, see the full proof in the following section.To prove Theorem 1(b), it suffices to show that the moment-generating function (MGF) of the QAOA overlap converges to the MGF of a sine-Gaussian law as follows:lim n→∞ E Y [M n (ζ)] = E G∼N (0,1) exp ζe -2qγ 2 sin(2β) sin(2qΛγG q-1 ) =: M (ζ). (C.4)We start the proof of Eq. (C.4) with the following lemma, which obtains a more explicit expression for the MGF that we derived in Section B.Lemma C.1 (Expected moment-generating function). The expected moment-generating function in Eq. (B.4) can be evaluated asE Y [M n (ζ)] = n t=0 n t e -γ 2 [n q -(n-2t) q ]/n q-1 [sinh(ζ/n) sin(2β)] t • E n,t , (C.5)whereE n,t = 1 2t + 1 t ξ=-t sin t (2πξ/(2t + 1)) Ẑn,t (ξ), Ẑn,t (ξ) = t k=-t e -2πiξk/(2t+1) Z n,t (k), Z n,t (k) = 1 2 n-t τ++τ-=n-t n -t τ + , τ - (e ζ/n cos 2 β + e -ζ/n sin 2 β) τ+ (e -ζ/n cos 2 β + e ζ/n sin 2 β) τ-× e iΛnγ[((τ+-τ-)+k) q -((τ+-τ-)-k) q ]/n (q-1)/2 . (C.6) Here we haveΛ n = λ n /n (q-1)/2 .The proof of Lemma C.1 is deferred to Section C.2.1. Now we defineΛ = lim n→∞ Λ n , I n,t = n t e -γ 2 [n q -(n-2t) q ]/n q-1 [sinh(ζ/n) sin(2β)] t • E n,t , I t = 1 t! E G∼N (0,1) [ζe -2qγ 2 sin(2β) sin(2qΛγG q-1 )] t . (C.7)Then it is easy to see thatE Y [M n (ζ)] = n t=0 I n,t , M (ζ) = ∞ t=0 I t .As a consequence, we haveE Y [M n (ζ)] -M (ζ) ≤ T t=0 |I n,t -I t | + t≥T +1 I t + n t=T +1 |I n,t |. (C.8)The following lemma gives the limit of E n,t for fixed t as n → ∞, which indicates that I t is the limit of I n,t . Lemma C.2. For any fixed integer t, we havelim n→∞ E n,t = E G∼N (0,1) [sin t (2qΛγG q-1 )] ≡ E t .As E Y [M n (ζ)] -M (ζ) ≤ Tε t=0 |I n,t -I t | + t≥Tε+1 I t + ∞ t=Tε+1 s t ≤ ε. (C.9)This proves Eq. (C.4) as desired, and hence finishes the proof of Theorem 1(b).Our starting point is Eq. (B.4), which we can compute explicitly with a careful organization of the sum. To this end, lett + = n ++-+ n -++ , t -= n +--+ n --+ , d + = n ++--n -++ , d -= n +---n --+ , τ + = n +++ + n ---, τ -= n +-+ + n -+-, ∆ + = n +++ -n ---, ∆ -= n +-+ -n -+-.(C.10)Observe that these 8 variables completely determine {n a : a ∈ B}. Furthermore, lett = t + + t -, n -t = τ + + τ -. (C.11)Then explicit computation shows thata∈B q Φ 2 a q s=1 n as = 4γ 2 a 1{a 11 • • • a q1 ̸ = a 12 • • • a q2 } q s=1n aq = 2γ 2 n q -(n -2t) q , (C.12)a∈B q Φ a q s=1 (a s ) m n as = γ a a 1 a m n a q - a a 2 a m n a q = γ[(τ + -τ -) + (d + -d -)) q -((τ + -τ -) -(d + -d -)) q ], (C.13) v∈B v m n v = t + -t -+ ∆ + -∆ -. (C.14)Plugging this into Eq. (B.4) and breaking up the sum, we getE Y [M n (ζ)] = n t=0 n t e -γ 2 [n q -(n-2t) q ]/n q-1 t++t-=t t t + , t -τ++τ-=n-t n -t τ + , τ - ∆+ τ + n +++ Q n+++ +++ Q n--- --- ∆- τ - n +-+ Q n+-+ +-+ Q n-+- -+-e (ζ/n)(t+-t-+∆+-∆-) d+ t + n ++- Q n++- ++-Q n-++ -++ d- t - n +-- Q n+-- +--Q n--+ --+ e iΛnγ[(d+-d-+τ+-τ-) q -((τ+-τ-)-(d+-d-)) q ]/n (q-1)/2 , (C.15)where Λ n = λ n /n (q-1)/2 as shorthand. Now, let us evaluate the sum over ∆ ± . We can use the following identity2 τ+ ∆+ τ + n +++ Q n+++ +++ Q n--- ---e (ζ/n)∆+ = (2Q +++ e ζ/n + 2Q ---e -ζ/n ) τ+ . (C.16)Applying this to the earlier sum, and note that Q +++ = Q +-+ and Q ---= Q -+-, we get:E Y [M n (ζ)] = n t=0 n t e -γ 2 [n q -(n-2t) q ]/n q-1 t++t-=t t t + , t -τ++τ-=n-t n -t τ + , τ - 1 2 n-t (2Q +++ e ζ/n + 2Q ---e -ζ/n ) τ+ (2Q +++ e -ζ/n + 2Q ---e ζ/n ) τ-e (ζ/n)(t+-t-) d+ t + n ++- Q n++- ++-Q n-++ -++ d- t - n +-- Q n+-- +--Q n--+ --+e iΛnγ[(d+-d-+τ+-τ-) q -((τ+-τ-)-(d+-d-)) q ]/n (q-1)/2 . (C.17)To further simplify the expression, we define Z n,t (k) as in Eq. (C.6), which we reproduce here:Z n,t (k) = 1 2 n-t τ++τ-=n-t n -t τ + , τ - (e ζ/n cos 2 β + e -ζ/n sin 2 β) τ+ (e -ζ/n cos 2 β + e ζ/n sin 2 β) τ- × e iΛnγ[((τ+-τ-)+k) q -((τ+-τ-)-k) q ]/n (q-1)/2 . (C.18)Then we haveE Y [M n (ζ)] = n t=0 n t e -γ 2 [n q -(n-2t) q ]/n q-1 t++t-=t t t + , t -d+ t + n ++- Q n++- ++-Q n-++ -++ × d- t - n +-- Q n+-- +--Q n--+ --+ e (ζ/n)(t+-t-) Z n,t (d + -d -).(C. 19 ) Let Ω t = {-t, -t + 1, . . . , t -1, t}, and let { Ẑn,t (ξ)} ξ∈Ωt to be the discrete Fourier transform of Z n,t (k) as defined in Eq. (C.6), i.e.,Ẑn,t (ξ) = (F t Z n,t )(ξ) = t k=-t e -2πiξk/(2t+1) Z n,t (k), (C.20)By the property of Fourier transforms, we have Z n,t (k) = (F -1 t Ẑ)(k) = 1 2t + 1 t ξ=-t e 2πiξk/(EQUATION) (i) = n t=0 n t e -γ 2 [n q -(n-2t) q ]/n q-1 t++t-=t t t + , t - e (ζ/n)(t+-t-) × (-1) t-• 1 2t + 1 t ξ=-t 2iQ ++-sin(2πξ/(2t + 1)) t Ẑn,t (ξ) (C.23) (ii) = n t=0 n t e -γ 2 [n q -(n-2t) q ]/n q-1 (sinh(ζ/n) sin(2β)) t × 1 2t + 1 t ξ=-t sin(2πξ/(2t + 1))t Ẑn,t (ξ), (C. 24) where (i) used the equationd+ t + n ++- Q n++- ++-Q n-++ -++ e iξd+ = 2iQ ++-sin ξ t+ , (C.25)and (ii) used the equationr+s=t t r, s (+1) r (-1) s exp{ζ(r -s)} = 2 t sinh(ζ) t . (C.26)Note we have also used 4iQ ++-= sin 2β to get rid of the two factors of 2 t . This completes the proof of Lemma C.1.We first look at the limit of Z n,t (k) for fixed integer -t ≤ k ≤ t (c.f. Eq. (C.6)). We denoteT n = (e ζ/n cos 2 β + e -ζ/n sin 2 β), U n = (e -ζ/n cos 2 β + e ζ/n sin 2 β) and G n = (τ + -τ -)/ √ n with τ + ∼ Bin(n -t, 1/2)) to be a random variable. Then we haveZ n,t (k) = E Gn T √ nGn n (T n U n ) (n- √ nGn)/2 e iΛnγ √ n[(Gn+k/ √ n) q -(Gn-k/ √ n) q ]) . (C.27)Note that we have lim n→∞ T√ n n = lim n→∞ (T n U n ) - √ n/2 = lim n→∞ (T n U n ) n/2 =1 and by assumption we have lim n→∞ Λ n = Λ. Furthermore, by central limit theorem, we have G n → G ∼ N (0, 1) so that for any fixed-t ≤ k ≤ t, √ n[(G n + k/ √ n) q -(G n -k/ √ n) q ]) d -→ 2qkG q-1 .n→∞ Z n,t (k) = E G∼N (0,1) [e ik2qΛγG q-1 ] ≡ Z(k).As a consequence, we havelim n→∞ E n,t = 1 2t + 1 t ξ=-t sin(2πξ/(2t + 1)) t t k=-t e -2πiξk/(2t+1) E G∼N (0,1) [e ik2qΛγG q-1] .Finally, by Lemma C.4 below and noting that sin(2πξ/(2t + 1)) t can be expressed as a degree t polynomial of (e 2πiξ/(2t+1) , e -2πiξ/(2t+1) ), the right hand side of the equation above gives1 2t + 1 t ξ=-t sin(2πξ/(2t+1)) t t k=-t e -2πiξk/(2t+1) E G∼N (0,1) [e ik2qΛγG q-1 ] = E G∼N (0,1) [sin(2qΛγG q-1 ) t ].This proves Lemma C.2. Lemma C.4. Let t ∈ Z ≥0 be an integer and let Ω t = {-t, -t + 1, . . . , t -1, t}. For a vector (Z(k)) k∈Ωt , we denote F t : C 2t+1 → C 2t+1 to be the discrete Fourier transform(F t Z)(ξ) ≡ t k=-te -2πiξk/(2t+1) Z(k).Let P : C 2 → C be any fixed polynomials with degree less or equal to t ∈ Z ≥0 . Let X be a real-valued random variable. Then we have1 2t + 1 t ξ=-t P (e 2πiξ/(2t+1) , e -2πiξ/(2t+1) ) F t E X [e ikX ] (ξ) = E X [P (e iX , e -iX )]. (C.28)Proof of Lemma C.4. By linearity of the expectation operator and the discrete Fourier transform operator, we just need to prove Eq. (C.28) for P (e 2πiξ/(2t+1) , e -2πiξ/(2t+1) ) = e 2πipξ/(2t+1) for some integer -t ≤ p ≤ t. Note that we have1 2t + 1 t ξ=-t e 2πipξ/(2t+1) F t E X [e ikX ] = 1 2t + 1 t ξ=-t t k=-t e 2πipξ/(2t+1) e -2πikξ/(2t+1) E X [e ikX ] = E X 1 2t + 1 t ξ=-t t k=-t e 2πi(p-k)(ξ/(2t+1)-X/(2π)) e ipX = E X [e ipX ],(C. 29) where the last equality used the fact that1 2t + 1 t ξ=-t t k=-t e 2πi(p-k)(ξ/(2t+1)-X/(2π)) = 1for any integer -t ≤ p ≤ t and any real X. This completes the proof of Lemma C.4.By the definition of Z n,t (k) as in Eq. (C.6), it is easy to see that Proof of Lemma C.5. Denote Λ n = λ n /n (q-1)/2 . We can write|Z n,t (k)| ≤ 1 2 n-t τ++τ-=n-t n -t τ + , τ - e ζ/n cos 2 β + e -ζ/n sin 2 β τ+ e -ζ/n cos 2 β + e ζ/n sin 2 β τ- × e iΛnγ[((τ+-τ-)+k) q -((τ+-τ-)-k) q ]/n (q-1)/2 ≤ 1 2 n-t τ++τ-=n-t n -t τ + , τ - e τ+|ζ|/n e τ-|ζ|/n • 1 = e (n-E Y [⟨R 2 QAOA ⟩ γn,βn ] = ∂ 2 ∂ζ 2 ζ=0 E Y [M n (ζ; γ n , β n , λ n )]. (C.33)Using Eq. (C.5), we can see that only a few terms depend on ζ, whose derivative gives∂ 2 ζ ζ=0 sinh(ζ/n) sin(2β) t e ζ/n cos 2 (β) + e -ζ/n sin 2 (β) τ+ e -ζ/n cos 2 (β) + e ζ/n sin 2 (β) τ- = δ t=0 2n 2 2t sin(2β) + (τ + -τ -) 2 -(τ + + τ -) cos(4β) + (τ + -τ -) 2 + (τ + + τ -) + δ t=1 n 2 t(τ + -τ -) sin(4β) + δ t=2 n 2 t(t -1) sin 2 (2β). (C.34)Hence, only the t = 0, 1, 2 terms survive, and we can write 35 ) whereE Y [⟨R 2 QAOA ⟩ γn,βn ] = T 0 + T 1 + T 2 (C.T 0 = 1 2 n+1 n 2 τ++τ-=n n τ + , τ - (τ + -τ -) 2 -(τ + + τ -) cos(4β) + (τ + -τ -) 2 + (τ + + τ -) , (C.36) T 1 = sin(4β n ) 3n • 2 n-1 e -γ 2 n [n q -(n-2) q ]/n q-1 ξ∈{±1} sin(2πξ/3) 1 k=-1 e -2πiξk/3 × τ++τ-=n-1 n -1 τ + , τ - (τ + -τ -)e iΛnγn[((τ+-τ-)+k) q -((τ+-τ-)-k) q ]/n (q-1)/2 , (C.37) T 2 = (n -1) sin 2 (2β n ) 10n • 2 n-2 e -γ 2 n [n q -(n-4) q ]/n q-1 ξ∈{±1,±2}sin 2 (2πξ/5) n -2 τ + , τ - e iΛnγn[((τ+-τ-)+k) q -((τ+-τ-)-k) q ]/n (q-1)/2 . (C.38)Lemma C.5 then immediately follows from the Lemma C.7 below.Lemma C.7. For any n ≥ n 0 for some large n 0 , we haveT 0 = (1 + cos(4β n ))/(2n), |T 1 | ≤ 2 sin(4β n )e -qγ 2 n , |T 2 | ≤ 2 sin 2 (2β n )e -qγ 2 n .Proof of Lemma C.7. We can compute the first term directly as follows. Note thatτ++τ-=n n τ + , τ - (τ + -τ -) q = ∂ q ∂x q x=0 τ++τ-=n n τ + , τ - e x(τ+-τ-) (C.39) = ∂ q ∂x q x=0 2 cosh(x) n . (C.40)In particular,τ++τ-=n n τ + , τ - (τ + -τ -) = 0, (C.41) τ++τ-=n n τ + , τ - (τ + -τ -) 2 = 2 n n. (C.42)It follows thatT 0 = 1 2 n+1 n 2 ((2 n n -0) cos(4β) + (2 n n + 0)) = cos(4β) + 1 2n . (C.43)For the remaining terms, upper bounds suffice:|T 1 | ≤ sin(4β n ) 3n • 2 n-1 e -γ 2 n [n q -(n-2) q ]/n q-1 ξ∈{±1} 1 • 1 k=-1 1 • τ++τ-=n-1 n -1 τ + , τ - • (n -1) • 1 ≤ sin(4β n ) 3 e -qγ 2 n • 2 • 3 = 2 sin(4β n )e -qγ 2 n , (C.44)and D Derivation for general p-step QAOA (Claim 3.7)EQUATION• 2 k=-2 1 • τ++τ-=n-2 n -2 τ + , τ - • 1 ≤ sin 2 (2β n ) 10 e -qγ 2 n • 4 • 5 = 2 sin 2 (2β n )e -We now briefly sketch some ideas behind the derivation for Claim 3.7 that characterizes the overlap distribution of the p-step QAOA when the SNR ratio scales as in Eq. (3.10). Similar to Theorem 1, our approach is to evaluate the moment-generating function of the QAOA overlap in the n → ∞ limit. As evident in the proof of Theorem 1, as well as in previous analyses of the QAOA applied to spin-glass models [24, 26, 28] , the key technical difficulty is handling a "generalized multinomial sum" of the following form:EQUATIONwhere P (m) is a polynomial over entries of m = (m j ) j with degree q. Note the above summation has no analytical simplification when P is not a linear polynomial (q > 1). Previous works have evaluated this sum in the n → ∞ limit either by proving a "generalized multinomial theorem" that exploits combinatorial structures of the polynomial P [24, 26] , or by employing a Gaussian integration trick and the saddle-point method when q = 2 ℓ [28]. However, neither approach is sufficient for the spiked tensor model that we study in the present paper.Instead, we develop an alternative approach based on the Fourier transform to linearize exponents in the summands. The idea is to replace m with continuous variables µ via Dirac delta functions, which after Fourier transforms yield exponents that are linear in m, enabling us to analytically evaluate the multinomial sum over m as follows:S = ˆdµ ˆdμ mj ≥0, j mj =n n {m j } j Q mj j exp[P (µ)]e i μ•(m-µ) = ˆdµ ˆdμ j Q j e iμj n e P (µ)-i μ•µ . (D.2)See Appendix D.4 for more details. This is a powerful approach to replace the cumbersome multinomial sums with simpler integrals. However, it is difficult to make such manipulations involving Dirac delta functions rigorous, which we leave open as future work. Nevertheless, we proceed with the heuristic derivation in the current paper: by writing the variables (m j ) j in an alternative basis and rescaling them cleverly, we are able to evaluate the integrals to obtain the moment-generating function in the n → ∞ limit.Our goal is to evaluate the moment-generating function of the overlap with signal, M n (ζ) = ⟨γ, β| exp(ζ R)|γ, β⟩, for the general p-step QAOA. Using the same method as in the p = 1 case, we can show that the disorder-averaged moment-generating function can be written as the following combinatorial sum:E Y [M n (ζ)] = {na} n {n a } a∈B Q na a exp A + iλ n B + ζC , (D.3)whereA = - 1 2n q-1 a∈B q Φ 2 a q s=1 n as , B = 1 n q-1 a∈B q Φ a q s=1 (a s ) m n as , C = 1 n v∈B v m n v , (D.4)and B = (a 1 , a 2 , . . . , a p , a m , a -p , . . . , a -1 ) :a j ∈ {±1} , Q a = 1 2 p r=1 (cos β r ) 1+(ar+a-r)/2 (sin β r ) 1-(ar+a-r)/2 (i) (a-r-ar)/2 , Φ a = p r=1 γ r a r a r+1 • • • a p -a -p • • • a -r-1 a -r , Φ a = Φ a1a2•••aq . (D.5)Note Q a and Φ a are independent of a m . This is a straightforward generalization of the proof in Appendix B, where we insert 2p+1 resolutions of the identity instead of 3. This also closely follows the derivation in Ref. [26, Appendix D.2] . B, Q a , Φ a are also generalizations of the same quantities in Appendix B for p > 1.Define the rank function ℓ(a) = max({i :a -i ̸ = a i } ∪ {0}). (D.6)A canonical basis. We next perform further simplifications that remove the explicit dependence on a m . First we define the set of 2p-bit strings as A = (a 1 , a 2 , . . . , a p , a -p , . . . , a -1 ) : a j ∈ {±1} , and define A 0 and D according to a similar convention as that in [26] as follows:A 0 := {a ∈ A : ℓ(a) = 0} = {a ∈ A : a -k = a k for 1 ≤ k ≤ p}, D := a ∈ A : ℓ(a) > 0 and p j=1 a j = +1 . (D.7)Given the rank function in Eq. (D.6), we can define an ordering on D, which we borrow from [26] .For any two distinct element a 1 , a 2 ∈ D, we define the ≺ relation as following: (1) If ℓ(a 1 ) < ℓ(a 2 ), we let a 1 ≺ a 2 ; (2) If ℓ(a 1 ) > ℓ(a 2 ), we let a 2 ≺ a 1 ; (3) If ℓ(a 1 ) = ℓ(a 2 ) and if a 1 is lexically less than a 2 , we let a 1 ≺ a 2 ; (3) If ℓ(a 1 ) = ℓ(a 2 ) and if a 1 is lexically greater than a 2 , we let a 2 ≺ a 1 (here lexical order means that, for example, (-1, -1), (-1, 1), (1, -1), (1, 1) are in lexically increasing order). It is easy to see that such ≺ relation is a full order, so that we can also define ⪯, ⪰, and ≻ accordingly.For any a ∈ A, we define n a± = n b where b = (a 1 , . . . , a p , ±1, a -p , . . . , a -1 ). (D.8)Let We will call the last line as the "canonical basis". As a side note, comparing to the p = 1 derivation in Eq. (C.6), we have k = δd +-and τt a+ = n a+ + n ā+ , t a-= n a-+ n ā-, ∀a ∈ D, d a+ = n a+ -n ā+ , d a-= n a--n ā-, ∀a ∈ D, n a = n a+ + n a-, δn a = n a+ -n a-, ∀a ∈ A 0 . (D.+ -τ -= δn ++ -δn --.In what follows, we will convert all our expressions into the canonical basis. It is also helpful to denote the shorthand t = a∈D t a , and thus n -t = a∈A0 n a . (D.12)In this basis, we can rewrite (D.3) as 13) where we have used the fact that Q a± = Q a does not depend on a m (here we also slightly abused notation allowing Q a to take a ∈ A as argument). Here we also define, for any a ∈ D and t a ∈ Z ≥0 , the little-sum operator on functions of (d a , δd a , δt a ) as " ta da,δda,δtaE Y [M n (ζ)] = n t=0 n t {na} a∈A 0 n -t {n a } a∈A0 Q na a δna n a n a+ × {ta} a∈D t {t a } a∈D " ta da,δda,δta exp A + iλ n B + ζC . (D.(• • • ) := ta+,ta- t a t a+ , t a-da+ t a+ n a+ Q na+ a Q nā+ ā da- t a- n a- Q na- a Q nā- ā (• • • ). (D.14)Now let us rewrite B in the canonical basis, and we will show that it is purely a function of {δd a , δt a } a∈D ∪ {δn c } c∈A0 . Observe that To reveal additional structures of B, we writen q-1 B = p r=1 γ r B + r q -B -n q-1 B = p r=1 γ r [(R r + L r ) q -(R r -L r ) q ] = p r=1 γ r [2qL r R q-1 r + 2 q 3 L 3 r R q-3 r + • • • ] (D.16)where we have defined We note here that B consists of terms that have at least one power of the {δd a } a∈D variables through the dependence on L r , which is a fact that will become important later.L r = 1 2 (B + r -B - r ) = a∈D,ℓ(a)≥r 1 2 (a * r -a * -r )δd a , (D.17) R r = 1 2 (B + r + B - r ) =Proceeding in the same way for A and C, we can also write them in the canonical basis. We note A is a polynomial that has appeared in [26] , where it can be shown to only depend on {t a , d a } a∈D ∪ {n c } c∈A0 . In summary, we note the dependence of A, B, C on the canonical basis variables is as follows:A = A {t a } a∈D , {d a } a∈D , {n c } c∈A0 , iλ n B = iλ n B {δd a , δt a } a∈D ∪ {δn c } c,∈A0 , C = 1 n a∈A0 δn a + a∈D δt a . (D.19)Operator shorthands for different parts of the sum. To streamline notations, we now introduce three operators T, S, U as shorthands for different parts of the sum that appear in Eq. (D.13).Let us define the T t n operator acting on a function f ({t a : a ∈ D}) as ) and the summand A + iλ n B + ζC converge to simplified forms. To this end, for all a, b ∈ D, c ∈ A 0 , we will rescale by definingT t n f = t! n t n t ta≥0,t a = τ a , δt a /n ρa = δτ a , d b /n = η b , δd b /n 1-ρ b = δη b n c /n = ω c , δn c / √ n = δω c , (D.25)where (τ a , η b , ω c , δτ a , δη b , δω c ) are new dimensionless variables that will be integrated over, and ρ a are scaling exponents which we will define shortly.The goal of this subsection is to derive the summand in the n → ∞ limit. Specifically, we consider the summand broken into two parts, each as a polynomial of a distinct subset of the rescaled variables as follows:Γ n ({t a , η b , ω c }) := A({t a , η b n, ω c n}), (D.26) Ξ n ({δτ a , δη b , δω c }) := iλ n B({δτ a n ρa , δη b n 1-ρ b , δω c √ n}) + ζC({δτ a n ρa , δω c √ n}), (D.27)where the subscripts in the arguments implicitly iterate over a, b ∈ D and c ∈ A 0 . We think of Γ n and Ξ n as polynomials in their arguments, whose coefficients can depend on n.First, we know from [26, Lemma D.2] that with the rescaling specified in Eq. (D.25) and γ j , β j = Θ(1), we havelim n→∞ Γ n ({t a , η b , ω c } a,b∈D,c∈A0 ) = a∈D t a P a ({η b } b≺a , {ω c } c∈A0 ) =: Γ. (D.28)For the rest of this subsection, we derive the limit of Ξ n = iλ n B + ζC.Choosing the scaling exponents ρ a . We want to choose the scaling exponents for (δd a , δt a ) variables, such that all the terms of B except those are linear in δd a vanish in the n → ∞ limit. This would then imply the polynomial Ξ n in the limit would only be at most linear in δη a , which is very helpful later for evaluating certain integrals as we shall see in Eq. (D.61).In the general p-step QAOA applied to the spiked q-tensor model, suppose the SNR parameter λ has a scaling as follows λ n = Λn c(p,q) , (D. 29) where c(p, q) is to be determined. Also suppose that the appropriate scaling for δd a and δt a are a) , δt a ∼ n ρ ℓ(a) , (D.30) so that they only depend on the rank ℓ(a) of a. Based on the explicit derivation at p = 1, we believe we only care about the terms in B that look like δd a δn q-1 b when ℓ(a) = 1 and ℓ(b) = 0, or δd a δt q-1 b when ℓ(a) = ℓ and ℓ(b 25) . For these terms in B, we have λ n n q-1 δd a δn q-1 b ∼ n c(p,q)+1-ρ1+(q-1)/2-(q-1) , (D.31)δd a ∼ n 1-ρ ℓ() = ℓ -1 > 0. Also recall that δn b ∼ √ n for b ∈ A 0 from (D.λ n n q-1 δd a δt q-1 b ∼ n c(p,q)+1-ρ ℓ +(q-1)ρ ℓ-1 -(q-1) . (D.32)To ensure that all such terms in B are order 1, we impose the condition that c(p, q) + 1 -ρ ℓ + (q -1)(ρ ℓ-1 -1) = 0, andρ 0 = 1 2 . (D.33)Solving this recurrence equation, we get thatρ ℓ = 1 - (q -1) ℓ 2 + c(p, q) (q -1) ℓ -1 q -2 . (D.34)If we impose the additional condition that ρ p = 1 (so that δt a /n = Θ(1) to yield a nonvanishing overlap in C), this implies that the SNR scaling needs to bec(p, q) = q -2 2 (q -1) p (q -1) p -1 = q -2 2 + q -2 2[(q -1) p -1] . (D.35)Plugging this into Eq. (D.34), we getρ ℓ = 1 2 (q -1) p + (q -1) ℓ -2 (q -1) p -1 . (D.36)For the special case of q = 2, we have c(p, 2) = 1 2p , andρ ℓ = 1 2 + ℓ 2p .Note 1/2 ≤ ρ ℓ ≤ 1 since 1 ≤ (q -1) ℓ ≤ (q -1) p and 0 ≤ ℓ ≤ p. This means δd a = O(n 1/2 ) and δt a = Ω(n 1/2 ). Another property to note is that ρ ℓ is monotonically increasing with ℓ. In particular, ρ 0 = 1/2 and ρ p = 1.The limiting expression for Ξ n . To get the limiting polynomial for Ξ n = iλ n B + ζC, we substitute δd a = δη a n 1-ρa , δt a = δτ a n ρa , and δn c = δω c √ n, and take the n → ∞ limit. We first consider B as written in Eq. (D.16). In terms of the rescaled dimensionless variables, we haveL r = a∈D,ℓ(a)≥r 1 2 (a * r -a * -r )δη a n 1-ρa , R r = 1 2 (B + r + B - r ) = a∈A0 a * r δω a √ n + a∈D,ℓ(a)≤r-1 a * r δτ a n ρa + a∈D,ℓ(a)≥r 1 2 (a * r + a * -r )δη a n 1-ρa .With the exponents defined in Eq. (D.36), we note that L r is dominated by {δη a : ℓ(a) = r}, and R r is dominated by {δω a : a ∈ A 0 } when r = 1 and {δτ a : ℓ(a) = r -1} when r > 1. Thus, the appropriately rescaled L r and R r in the limit areEQUATIONRr := lim n→∞ R r n ρr-1 ≃ a∈A0 a * r δω a , r = 1 a∈D,ℓ(a)=r-1 a * r δτ a , r > 1 . (D.38)For λ n = Λn c(p,q) , we haveiλ n B = iλ n n q-1 p r=1 γ r k odd 2 q k L k r R q-k r , lim n→∞ iλ n B = lim n→∞ iΛn c(p,q) n q-1 p r=1 γ r k odd 2 q k Lk r Rq-k r n k(1-ρr)+(q-k)ρr-1 .One can verify that for any 1 ≤ r ≤ p,lim n→∞ n c(p,q) n q-1 n k(1-ρr)+(q-k)ρr-1 = n (k-1)(1-ρr-ρr-1) = 1, k = 1 1/n ϵ for some ϵ > 0, k ≥ 3 (D.39)Hence, in the n → ∞ limit, only the k = 1 term survives, andlim n→∞ iλB = iΛ p r=1 2qγ r Lr Rq-1 r . (D.40)Similarly, considerC = a∈A0 δn a n + a∈D δt a n = a∈A0 δω a √ n + a∈D δτ a n ρa n . (D.41)In the n → ∞ limit, the only terms that survive are δτ a when ℓ(a) = p for which ρ a = 1.Combining the two equations above, we have For succinctness, we denote the following vectors of (rescaled) variables t = (t a ) a∈D , d = (d b /n) b∈D , n = (n c /n) c∈A0 , δt = (δt a /n ρa ) a∈D , δd = (δd b /n 1-ρ b ) b∈D , δn = (δn c / √ n) c∈A0 . (D.43)We can then write the MGF asE Y [M n (ζ)] = E Y ⟨γ, β| exp(ζ 1 n n i=1 Z i )|γ, β⟩ = n t=0e n (t), (D. 44 )where e n (t) = T t n S {ta} n U t n exp Γ n (t, d, n) + Ξ n (δt, δd, δn) . (D.45)Here, Γ n and Ξ n are polynomials of their arguments whose coefficients can depend on n. Furthermore,T t n , S {ta} n, and U t n are summing operators defined in Eqs. (D.20), (D.21), (D.22) earlier. We now introduce dummy variables (δτ , η, δη, ω, δω) which will replace (δt, d, δd, n, δn) via Dirac delta functions:e n (t) = ˆδτ,η,δη,ω,δω T t n S {ta} n U t n exp Γ n (t, η, ω) + Ξ n (δτ , δη, δω) δ(δt -δτ )δ(d -η)δ(δd -δη)δ(n -ω)δ(δn -δω) = ˆδτ,η,δη,ω,δω ˆδτ,η,δ η, ω,δ ω T t n S {ta} n U t n exp Γ n (t, η, ω) + Ξ n (δτ , δη, δω) e iδ τ •(δt-δτ )+i η•(d-η)+iδ η•(δd-δη)+i ω•(n-ω)+iδ ω•(δn-δω) .where in the last line we used the Fourier representation of delta functions and introduced dual variables (δ τ , η, δ η, ω, δ ω).{ta} n is a sum over (d, δd, δt) and U t n is a sum over (n, δn). We can apply them directly to the relevant exponentials since their dependence is now linear, but involves the dual variables.First, let us evaluate the S {ta} n sum, which is defined in Eq. (D.21) as a composition of many littlesums. We start by considering a single little-sum with parameters (κ I , κ II , κ III ) of the following form:F a (κ I , κ II , κ III ) := " ta da,δda,δta e κIda+κIIδda+κIIIδta (D.46) = ta+,ta- t a t a+ , t a-da+ t a+ n a+ Q na+ a Q nā+ ā da- t a- n a- Q na- a Q nā- ā e κIda+κIIδda+κIIIδta .This can be evaluated using Q ā = -Q a and the basic identityda+ ta+ na+ (+1) na+ (-1) nā+ e κda+ = [2 sinh κ] ta+ . Applying this to the two inner sums in F a , we get thatF a (κ I , κ II , κ III ) = Q ta a ta+,ta- t a t a+ , t a- [2 sinh(κ I + κ II )] ta+ [2 sinh(κ I -κ II )] ta-e κIIIδta = (2Q a ) ta [sinh(κ I + κ II )e κIII + sinh(κ I -κ II )e -κIII ] ta = (4Q a ) ta (sinh κ I cosh κ II cosh κ III + cosh κ I sinh κ II sinh κ III ) ta . (D.47)Returning to S t n , we getS {ta} n [e iδ τ •δt+iδ η•δd+i η•d ] = a∈D n ta t a ! " ta da,δda,δta e iδ τ •δt+iδ η•δd+i η•d = a∈D (4nQ a ) ta t a ! i sin ηa n cos δ ηa n 1-ρa cos δτ a n ρa -cos ηa n sin δ ηa n 1-ρa sin δτ a n ρa ta . (D.48)Next, for U t n , we have from the multinomial theorem thatU t n [e i ω•n+iδ ω•δn ] = {na} a∈A 0 n -t {n a } a∈A0 Q na a δna n a n a+ e i ω•n+iδ ω•δn = a∈A0 2Q a e iωa/n cos δ ωa √ n n-t . (D.49)Take n → ∞ limit of e n (t). We now take the n → ∞ limit while keeping t fixed, assuming λ n = Λn c(p,q) . Recall the fact from Appendix D.3 that 0 < ρ a < 1 when ℓ(a) < p and ρ a = 1 when ℓ(a) = p. Then taking the n → ∞ limit of (D.48) yieldslim n→∞ S {ta} n [e iδ τ •δt+iδ η•δd+i η•d ] = a∈D (4Q a ) ta t a ! [g a (δτ a , ηa , δ ηa )] ta (D.50)whereg a (δτ a , ηa , δ ηa ) = iη a -δ ηa δτ a , ℓ(a) < p iη a cos δ ηa -δτ a sin δ ηa , ℓ(a) = p . (D.51)Similarly, taking the n → ∞ limit of (D.49) giveslim n→∞ U t n [e i ω•n+iδ ω•δn ] = exp a∈A0 2Q a (iω a - 1 2 δ ω2 a ) , (D.52)where we used the fact that a∈A0 2Q a = 1. We also note that for any sequence of functions {f n (t)} n that pointwise converges to f (t), we havelim n→∞ T t n f n (t) = lim n→∞ t! n t n t ta≥0,∀a∈D, a ta=t f n (t) = ta≥0,∀a∈D, a ta=t f (t) =: T t f (t). (D.53)Plugging these back into e n (t), we get in the limit e(t) := lim n→∞ e n (t) = ˆδτ,η,δη,ω,δω ˆδτ,η,δ η, ω,δ ω T t e Γ(t,η,ω)+Ξ(δτ ,δη,δω) e -iδ τ•δτ -i η•η-iδ η•δη-i ω•ω-iδ ω•δω e i ω•(2Q)-1 2 δ ω•(2Q δ ω) a∈D (4Q a ) ta t a ! [g a (δτ a , ηa , δ ηa )] ta .(D.54) where we denoted the vector Q = (Q a ) a∈A0 , and (2Q δω) j = 2Q j δω j to mean element-wise product.Sum over e(t) to get MGF. Now we perform the sum over t to get the moment-generating function of the overlap distribution, since (heuristically) lim n→∞ E Y [M n (ζ)] = ∞ t=0 e(t). Note that ∞ t=0 T t f (t) = ta≥0,a∈D f (t). (D.55)So in the n → ∞ limit, effectively we are summing over {t a } independently. We can also use the fact from [26, Lemma D.2] that Γ(t, η, ω) is linear in t,Γ(t, η, ω) = a∈D t a P a (η, ω). (D.56) Hence, we have ∞ t=0 e(t) = ˆδτ,η,δη,ω,δω ˆδτ,η,δ η, ω,δ ω e -iδ τ •δτ -i η•η-iδ η•δη e i ω•(2Q-ω) e -iδ ω•δω-1 2 δ ω•(2Q δ ω)exp a∈D 4Q a g a (δ τ , η, δ η)e Pa(η,ω) e Ξ(δτ ,δη,δω) .(D.57)The integrals over (ω, ω) yield Dirac delta functions that set each ω a = 2Q a . The integral over δ ω yields a Gaussian density function for δω, each with mean 0 and variance 2Q a . So we can set δω a = G a ∼ N (0, 2Q a ), and replace the integrals over (δ ω, δω) with an expectation over G = (G a ) a∈A0 . Our expression then simplifies to∞ t=0 e(t) = E G ˆδτ,η,δη ˆδτ,η,δ η e -iδ τ •δτ -i η•η-iδ η•δη exp a∈D 4Q a g a (δ τ , η, δ η)e Pa(η,2Q) e Ξ(δτ ,δη,G)=: E G ˆδτ,η,δη ˆδτ,η,δ η e S . (D.58)To do the remaining integrals, it is necessary to use additional structure of the polynomials P a , g a and Ξ. From [26] , we know there is an ordering (≺) of the elements of D such that the η dependence in P a is only on {η b : b ≺ a}. Furthermore, from Appendix D.3, we know Ξ(δτ , δη, δω) has a particular form:Ξ(δτ , δη, δω) = i a∈D δη a R a (δτ , δω) + ζ b:ℓ(b)=p δτ b . (D.59)We also know that the δτ dependence in R a is only on {δτ b : ℓ(b) < ℓ(a)}. More explicitly, from Eq. (D.38),R a (δτ , G) = 2qΛγ r a * r X q-1 r, where r = ℓ(a) andX r = b∈A0 b * r G b , r = 1 b∈D,ℓ(b)=r-1 b * r δτ b , r > 1 . (D.60)Let us now write out the exponent S in (D.58) using the form of g a in (D.51) and Ξ in (D.59): Integrating over (η, η) yields delta functions that assign η a = 4Q a e Pa(η,2Q) when ℓ(a) < p, or η a = 4Q a cos δ ηa e Pa(η,2Q) when ℓ(a) = p. Note 4Q a e Pa = 2W a where W a is defined the same way for q-spin models as in [26] , so we will use W a = 2Q a e Pa in what follows. Then, integrating over (δη, δ η) yields delta functions that assign δ ηa = R a (δτ , G). Note here the linear dependence in δη a in Ξ(δτ , δη, δω), as in Eq. (D.59), is important for allowing us to evaluate the integrals. Finally, integrating over (δ τ , δτ ) yields delta functions that assign δτ a = i4Q a R a e Pa = i2W a R a when ℓ(a) < p, and δτ a = i4Q a sin R a e Pa = i2W a sin R a when ℓ(a) = p. Note that these assignments by delta functions are consistent if we perform the integrals according to the ascending order of the set D, since P a , R a only depend on the variables {(η b , δτ b ) : b ≺ a}, which would have already been assigned values from earlier integrals.S = -iδ τ • δτ -iη • η -iδ η • δηThe MGF of the overlap distribution is thenlim n→∞ E Y [M n (ζ)] = ∞ t=0 e(t) = E G exp ζ b∈D:ℓ(b)=p i2W b sin R b (G) . (D.62)In what follows, let us denote D r = {b ∈ D : ℓ(b) = r} for 1 ≤ r ≤ p. Also let γr = 2qΛγ r , and G = a∈A0 a * 1 G a . Note G ∼ N (0, 1) since G a ∼ N (0, 2Q a ) and a∈A0 2Q a = 1. To get a sense of the MGF formula, observe thatℓ(a) = 1 =⇒ R a = γ1 a * 1 G q-1 , ℓ(a) = 2 =⇒ R a = γ2 a * 2 b∈D1 i2W b R b b * 2 q-1 = γ2 a * 2 b∈D1 i2W b b 1 q-1 [γ 1 G q-1 ] q-1 .Note in the last line we used b * 1 b * 2 = b 1 . Doing this iteratively, we see that when ℓ(a) = r, we haveR a = a * r K r G (q-1) r , where K r = γr b∈Dr-1 i2W b b r-1 q-1 K q-1 r-1 (D.63)with initial condition K 1 = γ1 . Note that K r ∼ Λ [(q-1) r -1]/(q-2) when q > 2 and K r ∼ Λ r when q = 2. Furthermore, using the fact that sin(aX) = a sin X when a ∈ {±1}, we have from Eq. (D.62) thatR QAOA d -→ a∈Dp i2W a a * p sin K p G (q-1) p , (D.64)which is indeed of the form of the sine-Gaussian law in Claim 3.7.We then note that the factors b∈Dr-12W b b r-1 , b∈Dp 2W b b * p (D.65)can be evaluated efficiently using the iterative procedure in [25] due to Theorem 3 in [26] . We give this procedure in the section that immediately follows. This concludes the derivation that shows Claim 3.7.D.5 A self contained formula for a p (γ, β) and b p (γ, β)In this section, we give a self-contained description of the formula for (a p , b p ), following Eq. (D.64). Let B be the set of (2p + 1)-bit strings indexed as B = (z 1 , z 2 , . . . , z p , z 0 , z -p , . . . , z -1 ) : z j ∈ {±1} . Definef (z) = 1 2 ⟨z 1 |e iβ1X |z 2 ⟩ • • • ⟨z p-1 |e iβp-1X |z p ⟩ ⟨z p |e iβpX |z 0 ⟩ × ⟨z 0 |e -iβpX |z -p ⟩ ⟨z -p |e -iβp-1X |z -(p-1) ⟩ • • • ⟨z -2 |e -iβ1X |z -1 ⟩ (D.66)where z i ∈ {+1, -1}, and ⟨z 1 |e iβX |z 2 ⟩ = cos β if z 1 = z 2 , or i sin(β) otherwise. Define matrices H [m] ∈ C (2p+1)×(2p+1) for 0 ≤ m ≤ p as follows. For j, k ∈ {1, . . . , p, 0, -p, . . . , -1}, letEQUATIONH [m] j,k = z∈B f (z)z j z k exp - q 2 p j ′ ,k ′ =-p H [m-1] j ′ ,k ′ q-1 γ j ′ γ k ′ z j ′ z k ′ for 1 ≤ m ≤ p, (D.67)where we use the convention that γ -r = -γ r for 1 ≤ r ≤ p, and γ 0 = 0. Note these matrices first appeared in [25] in the context of assessing the performance of the QAOA on locally treelike Max-q-XORSAT problems and can be evaluated in O(p 2 4 p ) time.Once we have the matrix H [p] , we compute for 1 ≤ r ≤ p, Example formula at p = 2. As an example, we now describe the explicit formula at p = 2, which applies in the regime where λ n = Λn (q-2+1/q)/2 (note here ε p=2 = 1/q). We have b 2 = 2 q q q-1 e -2q(q-1)γ 2 1 γ q-1 1 γ 2 sin q-1 (2β 1 ),a r = i z∈B f (z) z r z r+1 -z -r z -(r+1) 2 p s=r+1 1 + z s z -s 2 exp - q 2 p j,k=-pa 2 = -e -2q(γ 2 1 +γ 2 2 +2 Re[X]γ1γ2) sin 2β 2 × cos 2 β 1 + e 8qγ1γ2 Re[X] sin 2 β 1 + e 2q(γ 2 1 +2γ1γ2 Re[X]) sin 2β 1 sin(4qγ 1 γ 2 Im[X]) ,where X = (cos 2β 1 + ie -2qγ 2 1 sin 2β 1 ) q-1 . Then the overlap R d -→ a 2 sin(b 2 Λ q G (q-1) 2 ).Although the above formula is complicated, we can understand the scaling with q by considering a simple choice of γ 1 = γ 2 = 1/2 √ q and β 1 = β 2 = π/4. Then the above simplifies to b 2 = e (1-q)/2 q q/2 , a 2 = e -1 cosh[e (1-q)/2 sin(πq/2)] -e -1/2 sin[e (1-q)/2 cos(πq/2)].(D.70)E Signal boosting with 1-step QAOA Consider a scenario where we have some prior information about the signal, in the form of a weak estimator that overlaps partially with the true signal. Our goal is to boost the overlap of this weak estimator. We study the SNR threshold of the 1-step of QAOA and compare it to the 1-step of power iteration. For QAOA, we encode the weak estimator into the initial state: rather than initializing with the uniform superposition across all bit-strings |s⟩, we bias a fraction of the qubits toward the signal. For power iteration, instead of starting from a uniform vector, we sample from a Bernoulli distribution biased toward the signal.More precisely, for QAOA we consider the following initial state:|s biased ⟩ = n j=1 cos θ j |u j ⟩ + sin θ j |-u j ⟩ , (E.1)where the θ j are drawn i.i.d. according toθ j = π/4, with probability 1 -k n , π/4 -δ, with probability k n , (E.2)and δ > 0. As in Eq. (2.3), we prepare the 1-step QAOA state as |γ, β⟩ biased = e -iβB e -iγC |s biased ⟩.Note the spiked tensor model Y is encoded in this state through C(σ) = ⟨Y , σ ⊗q ⟩/n (q-2)/2 . The following theorem concerns the SNR threshold for weak recovery and the distribution of overlap R QAOA,biased = û⊤ u/n between a sample û ∼ |γ, β⟩ biased and the signal u.Theorem 2 (Signal boosting with 1-step QAOA). Consider the biased 1-step QAOA state |γ, β⟩ biased as defined above. Fix γ > 0, β ∈ [0, 2π], δ ∈ [0, π/4], and let k = Θ(n c ) for 1/2 < c < 1. Supposelim n→∞ λ n /n (1-c)(q-1) = Λ. (E.3)Then, over the randomness of θ, Y and quantum measurement, the overlap R QAOA,biased of 1-step QAOA converges in probability to R QAOA,biased p -→ e -2qγ 2 sin(2β) sin(2qΛγ sin q-1 (2δ)).(E.4)We give the proof Theorem 2 in Appendix E.1 that follows.Remark E.1. Theorem 2 considers an initial state with a fraction k/n of qubits biased toward the signal vector u, representing some side information. It shows that the SNR threshold is Θ(n (1-c)(q-1) ), which becomes lower with increasing side information k/n = n c-1 . In particular, if k = Θ(n 3/4 ), the weak recovery threshold of 1-step QAOA improves to Θ(n (q-1)/4 ), compared to the Θ(n (q-1)/2 ) threshold given by Theorem 1 without any initial overlap between the state and planted signal.Comparison with classical tensor power iteration. We compare the boosting produced by the 1-step QAOA to that provided by 1-step power iteration. Recall the 1-step tensor power iteration estimator (2.1) is û1,biased = √ nY û⊗(q-1) 0,biased / Y û⊗(q-1) 0,biased 2 , where in this case, analogously to Eq. (E.1), the initial vector û0,biased has its entry (û 0,biased ) j sampled as (û 0,biased ) j ∼ u j / √ n, with probability 1 2 1 + k n sin(2δ) , -u j / √ n, with probability 1 2 1 -k n sin(2δ) .(E.5)One can check that √ nû 0,biased ∼ |s biased ⟩ is a sample from the biased initial QAOA state, so that we are making a fair comparison with QAOA. In the following proposition, we show that the required SNR for the 1-step power iteration estimator is also Θ(n (1-c)(q-1) ), and we provide the distribution of overlap R PI,biased ≡ u ⊤ û1,biased /n between the power iteration estimator û1,biased and the signal u. Proposition E.2 (Signal boosting with 1-step tensor power iteration). Assume that the rescaled signal-to-noise ratio has a limit lim n→∞ λ n /n (1-c)(q-1) = Λ. Then over the randomness of W and initialization û0,biased , the overlap R PI,biased of the power iteration estimator with the signal converges in probability to R PI,biased p -→ sin[arctan(Λ sin q-1 (2δ))].(E.6)The proof of Proposition E.2 is contained in Appendix H.2. This shows yet again that the QAOA has the same asymptotic computational efficiency as power iteration. Nevertheless, in the Λ ≪ 1 regime, by choosing γ = 1/2 √ q and β = π/8, the QAOA achieves an overlap that is larger than power iteration by a factor q/e.Without loss of generality, we assume that u = 1. Recall that the initial state is given by Eq. (E.1), which we can rewrite asEQUATIONwhere θ j = π/4 with probability 1 -k/n, and θ j = π/4 -δ with probability k/n.To prove Theorem 2, it suffices to show that the moment-generating function (MGF) of the QAOA overlap converges to the MGF of a deterministic variable as follows:lim n→∞ E θ E Y [M n (ζ)] = exp ζe -2qγ 2 sin(2β) sin(2qΛγ sin(2δ) q-1 ) =: M (ζ). (E.8)The argument for the proof is the same as that for Theorem 1(b), except that we must prove analogous versions of Lemma C. E θ E Y [M n (ζ)] = n t=0 n t e -γ 2 [n q -(n-2t) q ]/n q-1 sinh(ζ/n) sin(2β) 1 - k n + k cos(2δ) n t • E n,t , (E.9) where E n,t = 1 2t + 1 t ξ=-t sin t (2πξ/(2t + 1)) Ẑn,t (ξ), Ẑn,t (ξ) = t l=-t e -2πiξl/(2t+1) Z n,t (l), Z n,t (l) = 1 2 n-t τ++τ-=n-t n -t τ + , τ - (e ζ/n cos 2 β + e -ζ/n sin 2 β) τ+ (e -ζ/n cos 2 β + e ζ/n sin 2 β) τ- × 1 + k sin(2δ) n τ+ 1 - k sin(2δ) n τ-× e iΛnγ[((τ+-τ-)+l) q -((τ+-τ-)-l) q ]/n c(q-1) . (E.10)The proof of Lemma E.3 is deferred to Section E.2. Note the only difference from the unbiased case (Lemma C.1) is the presence of the two terms1 - k n + k cos(2δ) n t and 1 + k sin(2δ) n τ+ 1 - k sin(2δ) n τ-, and the rescaled power of n in the exponent.We further defineΛ = lim n→∞ Λ n , I n,t = n t e -γ 2 [n q -(n-2t) q ]/n q-1 sinh(ζ/n) sin(2β) 1 - k n + k cos(2δ) n t • E n,t , I t = 1 t! [ζe -2qγ 2 sin(2β) sin(2qΛγ sin q-1 (2δ))] t , (E.11)where the definition of E n,t is given in Eq. (E.10). Then it is easy to see thatE θ E Y [M n (ζ)] = n t=0 I n,t , M (ζ) = ∞ t=0 I t .As a consequence, we haveE Y [M n (ζ)] -M (ζ) ≤ T t=0 |I n,t -I t | + t≥T +1 I t + n t=T +1 |I n,t |. (E.12)The following lemma gives the limit of E n,t for fixed t as n → ∞, which indicates that I t is the limit of I n,t . Lemma E.4. For any fixed integer t, we have The proof of Lemma E.4 and E.5 is deferred to Section E.3 and E.4, respectively. Now we assume that these two lemmas hold. By the fact that ∞ t=0 I t is finite and by Lemma E.5, for any ε > 0, there exists T = T ε such thatlim n→∞ E n,t = sin t (2qΛγ sin q-1 (2δ)) ≡ E t . (E.t≥Tε+1 I t ≤ ε/3, t≥Tε+1 s t ≤ ε/3.Furthermore, by Lemma E.4, there exists N = N ε such that as long as n ≥ N ε , we have As a consequence, by Eq. (E.12), for any n ≥ n ε and ζ ≤ n, we haveE Y [M n (ζ)] -M (ζ) ≤ Tε t=0 |I n,t -I t | + t≥Tε+1 I t + ∞ t=Tε+1 s t ≤ ε. (E.16)This proves Eq. (E.8) as desired, and hence finishes the proof of Theorem 2.With an added expectation over θ, Eq. (B.4) still holds with a modified Q a :E θ E Y [M n (ζ)] = {na} n {n a } a∈B Q na a exp - 1 2n q-1 a∈B q Φ 2 a q s=1 n as + iλ n n q-1 a∈B q Φ a q s=1 (a s ) m n as + ζ n v∈B v m n v , (E.17)andQ (a1,am,a2) = f β,k,δ (a 1 , a m , a 2 ) (E.18)with f defined below:f β,k,δ (z 1 j , z m j , z 2 j ) =                            sin 2 β, if (z 1 j , z m j , z 2 j ) = (-1, -1, -1). (E.19)This proof follows very closely that of Theorem 1(b) in Appendix C. From the change of variables in Eq. (C.10) to the breaking up in Eq. (C.15), the same expression still hold, except that we redefine Λ n = λ n /n (1-c)(q-1) , which amounts to the power of n in the exponential changing: when compared to Eq. (C.15):E θ E Y [M n (ζ)] = n t=0 n t e -γ 2 [n q -(n-2t) q ]/n q-1 t++t-=t t t + , t -τ++τ-=n-t n -t τ + , τ - ∆+ τ + n +++ Q n+++ +++ Q n--- --- ∆- τ - n +-+ Q n+-+ +-+ Q n-+- -+-e (ζ/n)(t+-t-+∆+-∆-) d+ t + n ++- Q n++- ++-Q n-++ -++ d- t - n +-- Q n+-- +--Q n--+ --+ e iΛnγ[(d+-d-+τ+-τ-) q -((τ+-τ-)-(d+-d-)) q ]/n c(q-1) . (E.20)However, it is not true anymore that Q +++ = Q +-+ and Q ---= Q -+-in general. We use the identity in Eq. (C.16) to writeE θ E Y [M n (ζ)] = n t=0 n t e -γ 2 [n q -(n-2t) q ]/n q-1 t++t-=t t t + , t -τ++τ-=n-t n -t τ + , τ - 1 2 n-t (2Q +++ e ζ/n + 2Q ---e -ζ/n ) τ+ (2Q +-+ e -ζ/n + 2Q -+-e ζ/n ) τ-e (ζ/n)(t+-t-) d+ t + n ++- Q n++- ++-Q n-++ -++ d- t - n +-- Q n+-- +--Q n--+ --+e iΛnγ[(d+-d-+τ+-τ-) q -((τ+-τ-)-(d+-d-)) q ]/n c(q-1) .Then we redefineZ n,t (l) = 1 2 n-t τ++τ-=n-t n -t τ + , τ - (e ζ/n cos 2 β + e -ζ/n sin 2 β) τ+ (e -ζ/n cos 2 β + e ζ/n sin 2 β) τ- 1 + k sin(2δ) n τ+ 1 - k sin(2δ) n τ-e iΛnγ[((τ+-τ-)+l) q -((τ+-τ-)-l) q ]/n c(q-1) (E.21)and, analogously to Eq. (C. 19 ) writeE θ E Y [M n (ζ)] = n t=0 n t e -γ 2 [n q -(n-2t) q ]/n q-1 t++t-=t t t + , t -d+ t + n ++- Q n++- ++-Q n-++ -++ × d- t - n +-- Q n+-- +--Q n--+ --+ e (ζ/n)(t+-t-) Z n,t (d + -d -).(E.22) Using the discrete Fourier transform, we haveE θ E Y [M n (ζ)] = n t=0 n t e -γ 2 [n q -(n-2t) q ]/n q-1 t++t-=t t t + , t -d+ t + n ++- Q n++- ++-Q n-++ -++ × d- t - n +-- Q n+-- +--Q n--+ --+ e (ζ/n)(t+-t-) 1 2t + 1 t ξ=-t e 2πiξ(d+-d-)/(2t+1) Ẑn,t (ξ) = n t=0 n t e -γ 2 [n q -(n-2t) q ]/n q-1 t++t-=t t t + , t - e (ζ/n)(t+-t-) × (-1) t-• 1 2t + 1 t ξ=-t 2iQ ++-sin(2πξ/(2t + 1)) t Ẑn,t (ξ) (E.23) since the same relations between Q ++-, Q -++ , Q +--, Q --+ hold. Finally, E θ E Y [M n (ζ)] = n t=0 n t e -γ 2 [n q -(n-2t) q ]/n q-1 (sinh(ζ/n) sin(2β)) t 1 - k n + k cos(2δ) n t × 1 2t + 1 t ξ=-t sin(2πξ/(2t + 1)) t Ẑn,t (ξ), (E.24)which is analogous to Eq. (C.24). This completes the proof of Lemma E.3.We first look at the limit of Z n,t (l) for fixed integer -t ≤ l ≤ t. Letting T n = (e ζ/n cos 2 β + e -ζ/n sin 2 β), U n = (e -ζ/n cos 2 β + e ζ/n sin 2 β) and ϵ = k sin(2δ)/n, we can writeZ n,t (l) = 1 2 n-t τ++τ-=n-t n -t τ + , τ - T τ+ n U τ- n (1 + ϵ) τ+ (1 -ϵ)τ-e iΛnγ[((τ+-τ-)+l) q -((τ+-τ-)-l) q ]/n c(q-1) .We letEQUATIONwhere τ + ∼ Binom(n -t, (1 + ϵ)/2) so that G n → G ∼ N (0, 1) by the central limit theorem sinceE τ+ [τ + -τ -] = ϵn -t(ϵ + 1) and Var τ+ [τ + -τ -] = (n -t)(1 -ϵ 2 ).Recall that ϵ = sin(2δ)n c-1 where 1/2 < c < 1. It follows that lim n→∞ T e -2πiξl/(2t+1) e iΛγ2l sin q-1 (2δ) = sin t (2qΛγ sin q-1 (2δ)), (E.29) where we used Lemma C.4 with X = 1 with probability 1. This completes the proof of Lemma E.4.((ϵ+1)n-t(ϵ+1))/2 n = lim n→∞ U (-(ϵ-1)n+t(ϵ+1))/2 n = 1 as well as lim n→∞ T √ n/2 n = lim n→∞ U √ n/2 n = 1. Hence, for any fixed -t ≤ l ≤ t, it follows that 1 n c(q-1) ( √ nG n + ϵn -t(ϵ + 1) + l) q -( √ nG n + ϵn -t(ϵ + 1) -l) q = 1 n c(q-1) ( √ nG n + n c sin(2δ) -t(n c-1 sin(2δ) + 1) + l) q -( √ nG n + n c sin(2δ) -t(n c-1 sin(2δ) + 1) -l) q → 2ql sin q-1 (2δWe first bound Z n,t (k):|Z n,t (l)| ≤ 1 2 n-t τ++τ-=n-t n -t τ + , τ - e ζ/n cos 2 β + e -ζ/n sin 2 β τ+ e -ζ/n cos 2 β + e ζ/n sin 2 β τ- × 1 - k sin(2δ) n τ+ 1 + k sin(2δ) n τ-e iΛnγ[((τ+-τ-)+l) q -((τ+-τ-)-l) q ]/n c(q-1)≤ 1 2 n-t τ++τ-=n-t n -t τ + , τ - 1 - k sin(2δ) n τ+ 1 + k sin(2δ) n τ- e τ+|ζ|/n e τ-|ζ|/n • 1 = e (n-t)|ζ|/n • 1 ≤ e |ζ| .(E.30) The rest of the proof is exactly as the proof of Lemma C.3, except that I n,t involves the following extra factor which we can bound:1 - k n + k cos(2δ) n ≤ 3. (E.31)So the end bound on |I n,t | ends up with a different constant factor:|I n,t | ≤ 1 t! (18|ζ|) t (2t + 1)e |ζ| . (E.32)This finishes the proof of Lemma E.5.F Finite n calculation for 1-step QAOA on the spiked matrix (q = 2)In this appendix, we calculate the average squared overlap outputted by the QAOA at any finite problem dimension n and obtain the formula we reported in Eq. (4.1). As done in Appendix B, we first take u = 1 to be the all-one vector without loss of generality. The cost function isC(z) = n j,k=1 Y j,k z j z k , where Y j,k = λ n n + 1 √ n W j,k . (F.1)Here W j,k ∼ N (0, 1).The QAOA state at level p = 1 with this cost function is |γ, β⟩ = e -iβB e -iγC |s⟩ . (F.2)We are interested in the overlap of the QAOA output with the hidden signal u = 1. Following the same method as in Appendix B, we can write the disorder-averaged overlap as n t e -4γ 2 t(n-t)/n t++t-=t t t + , t -τ++τ-=n-t n -t τ + , τ -E Y [⟨R 2 QAOA ⟩ γ,β ] = {na} n {n a } a∈B Q na a e -1∆+ τ + n +++ Q n+++ +++ Q n--- --- ∆- τ - n +-+ Q n+-+ +-+ Q n-+- -+-(t + -t -+ ∆ + -∆ -) 2 d+ t + n ++- Q n++- ++-Q n-++ -++ d- t - n +-- Q n+-- +--Q n--+--+ e iΛ(d+-d-)(τ+-τ-) , (F.12)where we've denoted Λ = 4λγ/n as shorthand. Note we need to perform these sums in a carefully chosen order in order to get a closed-form answer at the end.We start with the last line, where we sum over d ± . We can use the fact that Q +x-= -Q -y+ = -i 2 sin β cos β for any x, y ∈ {±}. Then, for example we haved+ t + n ++- Q n++- ++-Q n-++-++ e iΛd+(τ+-τ-) = 2iQ Note the Kronecker deltas will collapse the sum over t, so it remains to evaluate the sums over ∆ ± and τ ± . To perform the sum over ∆ ± , note that Q +x+ = 1 2 cos 2 β and Q -y-= 1 2 sin 2 β, for any x, y ∈ {±}. Thus, we can use the following identity Thus, plugging in Λ = 4λγ/n and using the fact that 2iQ ++-= 1 2 sin 2β, we arrive at 2 τ+ ∆+ τ + n +++ Q n+++ +++ Q n--- ---(∆ + ) k =    1, if k = 0, τ + cos 2β, if k = 1, τ + [1 + (τ + -1) cos 2 (2β)], if k = 2, (F.E Y [⟨R 2 QAOA ⟩ γ,β ] = n -1 2n e -8γIn this appendix, complementing the simulation results in Section 4, additional numerical simulations are performed for 1 ≤ p ≤ 7 at q = 2 and 1 ≤ p ≤ 6 at q = 3.Fig. 4 displays the second moment of QAOA overlap versus problem dimension n. The y-axis plots the simulated second moment subtracting the theoretical value in the n → ∞ limit. For all demonstrated (p, q) pairs, the simulation appears to converge to the theoretical value with order 1/n deviations. This is consistent with the rigorous finite-n formula for (p, q) = (1, 2) in Eq. (F.24). Figure 4 : Log-log plots of the difference between observed overlap (averaged over instances and quantum measurements) at various problem dimension n and the predicted value from the sine-Gaussian law in the n → ∞ limit. Different colored lines correspond to different QAOA depth p, with parameters (γ, β) set to be the same as in Table 1 . We choose Λ = 0.2, λ n = Λn 1/(2p) (left), and λ n = Λn [1+1/(2 p -1)]/2 (right). Error bars are standard errors of the mean.Our numerical simulations use the GenQAOA library, available at https://github.com/ leologist/GenQAOA. The simulations are conducted on a laptop (MacBook Pro M2), where the simulation of each 26-qubit instance with p-step QAOA for 1 ≤ p ≤ 7 took about 160 seconds. Data used in the figures are available upon request.We take ε ∈ ([(q -1) p-1 -1]/[2(q -1) p -2], 1/2) to be fixed. By Lemma H.1, for a fixed p ∈ N + , with high probability, we have T ε ≥ p, as well as the upper bounds indicated in Eq. (H.5) for all t ≤ p -1. Applying Eq. (H.4) recursively implies that α p d = (1 + o P (1)) • n 1/2 Λ 1/εp G (q-1) p-1 , G ∼ N (0, 1).By the last equation on page 9 of [33] , we see that (for H p = {0, 1, 2, . . . , p} q-1 ) → sin[arctan(Λ 1/εp G (q-1) p )]. This concludes the proof of Proposition 3.9.v p = α p v +• At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.Another commonly studied prior is the uniform distribution over the n-sphere, Unif(S n-1 ( √ n)). In this work, we choose the Rademacher prior for the convenience of the QAOA analysis. 
 
 Appendix:  


Decision: Consider statistical estimation in the spiked tensor model. The *weak recovery threshold* refers to the signal-to-noise ratio (SNR) above which weak recovery—an estimator that achieves non-trivial overlap with the signal—is possible. The paper proposes computing the maximum-likelihood estimate via the quantum approximate optimization algorithm (QAOA) and analyzes the weak recovery threshold of $p$-step QAOA.

- It is rigorously proved that 1-step QAOA achieves the same weak recovery threshold as 1-step tensor power iteration.
- It is proved that, if a heuristically proved claim regarding the sine-Gaussian law of the overlap holds, then $p$-step QAOA can be a constant factor better than that of $p$-step tensor power iteration.

Although the results of this paper do not surpass state-of-the-art classical algorithms and are not fully rigorous, it marks the first theoretical work to study QAOA for statistical estimation in the spiked tensor model and provides useful insights into whether there is a quantum advantage for this task. Hence, I suggest accepting this paper as a spotlight.

Reviews: review: paper_summary: paper_summary: This paper studies the performance of 1-step and multi-step quantum approximate optimization algorithm (QAOA) for spiked tensor problem. In this problem one observes a q-dimentional tensor which is a properly normalized linear combination of q-th tensor power of unknown vector $u \in \{+1, -1\}^n$ and Gaussian noise $W$: $\lambda u^{\otimes q} / n^{q/2} + 1/sqrt(n)\cdot W$.  The goal is to recover the unknown vector $u$ from the observed tensor. This problem is known to be statistically solvable for $\lambda > T$, for some absolute constant T; however, it is known that under common complexity assumption, the problem has polynomial time classical algorithm only for $\lambda = \Omega(n^{(q-2)/4})$, providing a complexity gap.

This paper studies whether a specific family of quantum algorithms, called QAOA, can achieve a quantum advantage compared to classical algorithm. The paper proved negative result showing that constant-step QAOA applied to weak recovery problem in the spiked tensor model only achieves non-trivial overlap with the signal u when $\lambda = \Omega(n^{(q-1)/2})$ which nearly matches the threshold for the tensor power classical algorithm.
strength_weakness: strength_weakness: I think that the contribution of this paper is somewhat limited in the sense that negative result is obtained only for a very particular family of quantum algorithms, which does not rule out that even small modifications can achieve quantum advantage for this problem. This is in contrast to classical results where gap is established for any classical algorithm under some standard complexity assumptions.
questions: questions: 1) Do authors expect that under some standard assumptions, any BQP algorithm is not able to recover $u$ for $\lambda =o(n^{(q-2)/4})$?
2) Do you see other candidate problems where techniques developed in this paper can potentially be used to study the performance of QAOA?
limitations: limitations: na
paper_summary: paper_summary: This submission proposed to use quantum approximate optimization algorithm (QAOA) to compute the maximum likelihood estimator in the statistical estimation problem of the spiked tensor model.

Using the overlap between the estimated vector and the original vector, the author(s) obtained rigorous analysis for the so-called weak recovery threshold for 1-step QAOA. Namely, above such a threshold, the overlap will be non-zero with non-trivial probability; otherwise, the overlap will vanish with high probability. The author(s) also showed that the established weak recovery threshold matches the 1-step tensor power iteration classical algorithm.

For the p-step QAOA, the author(s) also obtained the weak recovery threshold based on some heuristic argument, which also matches that of the p-step tensor power iteration algorithm. Numerical experiments show that the QAOA method could achieve the state of the art threshold by combining with tensor unfolding (a technique used in the state of the art algorithm).

As far as I know, the submission proposed a new technical analysis for QAOA (by showing that the overlap exhibits an asymptotic sine-Gaussian distribution), providing a rigorous study of the polynomial-time QAOA. 
As such, I believe the work is worth sharing to the community. Hence, I recommend accept.

(I did not have enough time to check all the derivations provided in appendices in details, but the overall proof ideal seems logical to me.)
strength_weakness: strength_weakness: 1. The asymptotic analysis requires the number of qubits n approaching infinity, which is practically demanding.

2. The analysis of the p-step QAOA is based on an heuristic calculation.

3. It is not known if QAOA could achieve the same threshold as in the state of the art classical algorithm. 
Numerical experiments do provide some potentials of matching the result. Whether there is a quantum advantage remains unclear.
questions: questions: 1. In Remark 3.10, the author(s) claimed that for certain parameter p, the QAOA has a constant factor advantage over the classical power iteration algorithm in the overlap achieved.
The overlap via the p-step tensor is proved in Proposition 3.9. However, the overlap via QAOA given in (3.12) is based on heuristic calculations.
Hence, I'm not sure such an advantage is rigorous. If I did not misunderstand anything, the wording of claiming the constant factor advantage needs to be modified.

2. In Discussion, the author(s) asserted that "This implies that achieving a strong quantum advantage via the QAOA requires using a number of steps p that grows with n."
I understand that multiple (possibly infinite many) steps is needed to achieve a better performance for QAOA.
However, it is still not analytically evident if QAOA can match the state of the art threshold, because a heuristic calculation is required in the analysis.
Hence, whether there is a rigorous advantage for QAOA still remains unclear to me, even p approaching infinity.

3. As above, even if author(s) can rigorously show that p-step QAOA outperforms the p-step tensor power iteration, it probably not accurate to call it a "modest quantum advantage", since the p-step tensor power iteration is not the state of the art algorithm.
I highly recommend the author(s) to be more careful about the phrasing of quantum advantage throughout the paper.
limitations: limitations: The limitations are addressed in the manuscript. I also summarized them in Weaknesses.
Yet, since this submission is a theory work, I think the practical limitation is not a big concern.
paper_summary: paper_summary: The quantum approximate optimization algorithm is analyzed for the spiked tensor model. Weak recovery
of 1-step QAOA is rigorously shown to matche that of 1-step tensor power iteration. Heuristic calculations for p-step QAOA
matche that of p-step tensor power iteration.
strength_weakness: strength_weakness: For p-step QAOA the analysis is not rigorous. The observation of the intriguing sine-Gaussian law is numerical. Further analysis will be needed.
questions: questions: Maybe authors could discuss the limitations of implementing such algorithms on NISQ devices ? Any realistic prospects ?
limitations: limitations: Maybe authors could discuss the limitations of implementing such algorithms on NISQ devices ? Any realistic prospects ?
paper_summary: paper_summary: The paper investigates the performance of the Quantum Approximate Optimization Algorithm (QAOA) on the spiked tensor model problem. The authors demonstrate that QAOA's weak recovery threshold aligns with that of tensor power iteration and show through heuristic calculations that multi-step QAOA could potentially match but not exceed the classical computation threshold. A notable finding is the sine-Gaussian law for the asymptotic overlap distribution of p-step QAOA verified by simulations, which is distinct from classical methods and suggests a modest quantum advantage. The paper employs novel techniques, including Fourier transforms, to analyze the QAOA's performance and concludes with implications for potential quantum advantage in statistical inference problems.
strength_weakness: strength_weakness: - The results indicate that constant-step QAOA does not improve the recovery threshold beyond what is achievable by classical tensor power iteration by more than a constant factor, suggesting that the quantum advantage is modest.

- The paper does not address the performance of QAOA with more circuit depths, which is an open question and could be crucial for demonstrating a strong quantum advantage.
questions: questions: - Could the authors explain more about the generality of their proof techniques, i.e., could their proof techniques  be used in analysis for QAOA algorithms in other problem settings?
limitations: limitations: - The analysis for p-step QAOA (where p > 1) relies on heuristic arguments from physics, which may not be as rigorous as desired.
paper_summary: paper_summary: The paper studies the performance of the Quantum Approximate Optimization Algorithm (QAOA) for a classical average case problem from high dimensional statistic: tensor principal component analysis (tPCA), which exhibits a computational-statistical gap. The paper investigates if this algorithm can achieve a quantum advantage over its classical counterparts. The paper makes progress towards this question and suggests that the answer to expect is somewhat negative (only for QAOA). The main results are
1. After 1 step of QAOA, it archives weak recovery at the same SNR threshold (up to constants) as achieved by 1 step of clssical tensor power iteration. 
2. Using heuristic calculations (but not rigorously), they further showed that even after (some constant) $p$ steps of QAOA, the estimator succeeds at weak recovery at the same SNR as the tensor power iteration. This further suggests that even after using tensor unfolding, QAOA won't be able to surpass the computational threshold for this problem.
3. Along the way, they observe a sine Gaussian law for the asymptotic distribution of the overlap between the estimator (after $p$ steps of QAOA) and the ground truth. Again, this is proven for $p=1$ steps but empirically verified for $p>1$. 

Background: The problem is known to have a computational statistical gap. In particular, in the parameterization (1.1) taken in this paper, (a) recovery is possible whenever the SNR $\lambda \gg 1$, (b) but the threshold for efficient algorithms is known to be $\lambda \approx n^{(q-2)/4}$. Several classical algorithms including tensor unfolding, sum-of-squares, or gradient descent with landscape smoothing are known to achieve this. For iterative algorithms, such as the tensor power iteration and vanilla GD, the threshold is further away, requiring $\lambda \approx n^{(q-2)/2}$, and thus, to achieve the computational threshold either tensor unfolding (for power iteration) or landscape smoothing (for GD) is required.
strength_weakness: strength_weakness: I do not see any major weaknesses in the paper. Only a small quibble is a place in the introduction in lines 39-40 (and also in the abstract lines 1-3), where the authors present the motivation as seeking whether QAOA has superpolynomial speedup from clssical algorithms. However, I found this motivation slightly hand-wavy. I could not find enough concrete justification for why looking at (the combination of) QAOA for tensor PCA is a promising avenue for demonstrating this. If authors can make this more concrete, that would be helpful.

On the other hand, just studying the performance of QAOA for tensor PCA is an important question in its own right, as very well justified in lines 40-51. The authors do make good progress towards this.
questions: questions: 1. Could the author elaborate on how to combine QAOA with tensor unfolding? In more detail, why do the results after $p$ steps of unfolded tensors suggest that we could not surpass the computational threshold after using QAOA on unfolded almost square matrix?
2. In a non-quantum setup, a more standard is to take the prior to be uniform over a sphere. Can authors describe why it was important to take it to be uniform over the hypercube (in slightly more detail than in footnote 1)? I am just trying to understand the difficulty in the analysis if the prior was uniform over a sphere.
limitations: limitations: N/A

Comments: comment: We thank Reviewer a21g for their positive feedback and suggestion. In line with their recommendation, we will add more discussion about implementations of the QAOA in our revision.
We are grateful to Reviewer 9dK4 for their support. We will revise our manuscript according to their suggestions in the camera-ready version if this submission is accepted.
I thank the authors for their response, and I adjusted my score to 7.
Dear Reviewer TiV7, 

The author-reviewer discussion period is ending soon. Please check if the authors' response has addressed your concerns and feel free to adjust your score. If the authors' response is not satisfactory to you, please explain your reason and discuss with the authors *immediately*. 

Best regards, 
AC
Thanks for the detailed response! The authors have addressed my questions satisfactorily, so I adjusted the score accordingly.
Thanks for answering my questions! I would be happy to see this paper appear at the conference!
The authors have answered my question satisfactorily. I suggest adding a few pointers to the implementations of QAOA on curent devices to guide potentially interested readers. I understand this aspect in beyond the scope of the paper but if space allows a few pointers and comments woulb be welcome.
The authors have addressed my questions. Once the manuscript is revised accordingly, I can support this submission.