The Complexity of Sparse Tensor PCA
Year: 2021

Abstract:
We study the problem of sparse tensor principal component analysis: given a tensor Y = W + λx⊗p with W ∈ ⊗R having i.i.d. Gaussian entries, the goal is to recover the k-sparse unit vector x ∈ R. The model captures both sparse PCA (in its Wigner form) and tensor PCA. For the highly sparse regime of k ≤ √ n, we present a family of algorithms that smoothly interpolates between a simple polynomial-time algorithm and the exponential-time exhaustive search algorithm. For any 1 ≤ t ≤ k, our algorithms recovers the sparse vector for signal-tonoise ratio λ ≥ Õ( √ t · (k/t)) in time Õ(n), capturing the state-of-the-art guarantees for the matrix settings (in both the polynomial-time and sub-exponential time regimes). Our results naturally extend to the case of r distinct k-sparse signals with disjoint supports, with guarantees that are independent of the number of spikes. Even in the restricted case of sparse PCA, known algorithms only recover the sparse vectors for λ ≥ Õ(k ·r) while our algorithms require λ ≥ Õ(k). Finally, by analyzing the low-degree likelihood ratio, we complement these algorithmic results with rigorous evidence illustrating the trade-offs between signal-to-noise ratio and running time. This lower bound captures the known lower bounds for both sparse PCA and tensor PCA. In this general model, we observe a more intricate three-way trade-off between the number of samples n, the sparsity k, and the tensor power p.

text: √ n,
we present a family of algorithms that smoothly interpolates between a simple polynomial-time algorithm and the exponential-time exhaustive search algorithm. For any 1 ≤ t ≤ k, our algorithms recovers the sparse vector for signal-tonoise ratio λ ≥ Õ( √ t · (k/t)p/2) in time Õ(np+t), capturing the state-of-the-art guarantees for the matrix settings (in both the polynomial-time and sub-exponential time regimes). Our results naturally extend to the case of r distinct k-sparse signals with disjoint supports, with guarantees that are independent of the number of spikes. Even in the restricted case of sparse PCA, known algorithms only recover the sparse vectors for λ ≥ Õ(k ·r) while our algorithms require λ ≥ Õ(k). Finally, by analyzing the low-degree likelihood ratio, we complement these algorithmic results with rigorous evidence illustrating the trade-offs between signal-to-noise ratio and running time. This lower bound captures the known lower bounds for both sparse PCA and tensor PCA. In this general model, we observe a more intricate three-way trade-off between the number of samples n, the sparsity k, and the tensor power p.

heading: 1 Introduction
text: Sparse tensor principal component analysis is a statistical primitive generalizing both sparse PCA2 and tensor PCA3. We are given multi-linear measurements in the form of a tensor
Y = W + λx⊗p ∈ ⊗pRn (SSTM) for a Gaussian noise tensorW ∈ ⊗pRn containing i.i.d. N(0, 1) entries4 and signal-to-noise ratio λ > 0. Our goal is to estimate the “structured” unit vector x ∈ Rn. The structure we enforce on x is sparsity: |supp(x)| ≤ k. The model can be extended to include multiple spikes in a natural way: Y = W + ∑r q=1 λqx ⊗p (q), and even general order-p tensors: Y = W + ∑r q=1 λqX(q) for X(q) = x(q,1) ⊗ · · · ⊗ x(q,p) ∈ ⊗pRn. In this introduction, we focus on the simplest single spike setting of SSTM.
It is easy to see that sparse PCA corresponds to the setting with tensor order p = 2. On the other hand, tensor PCA is captured by effectively removing the sparsity constraint: |supp(x)| ≤ n. In ∗Part of the work was done while the author was in ETH Zürich. 2Often in the literature, the terms sparse PCA and spiked covariance model refer to the sparse spiked Wishart model. However, here we consider the sparse spiked Wigner matrix model. 3Tensor PCA is also known as the spiked Wigner tensor model, or simply the spiked tensor model. 4Throughout the paper, we will write random variables in boldface.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
recent years, two parallel lines of work focused respectively on sparse PCA [JL09, AW08, BR13a, DM16, HKP+17, DKWB19, HSV20, dKNS20] and tensor PCA [MR14, HSS15, MSS16, HKP+17, KWB19, AMMN19], however no result captures both settings. The appeal of the sparse spiked tensor model (henceforth SSTM) is that it allows one to study the computational and statistical aspects of these other fundamental statistical primitives in a unified framework, understanding the computational phenomena at play from a more general perspective.
In this work, we investigate SSTM from both algorithmic and computational hardness perspectives. Our algorithm improves over known tensor algorithms whenever the signal vector is highly sparse. We also present a lower bound against low-degree polynomials which extends the known lower bounds for both sparse PCA and tensor PCA, leading to a more intricate understanding of how all three parameters (n, k and p) interact.

heading: 1.1 Related work
text: Disregarding computational efficiency, it is easy to see that optimal statistical guarantees can be achieved with a simple exhaustive search (corresponding to the maximum likelihood estimator): find a k-sparse unit vector maximizing 〈Y , x⊗p〉. This algorithm returns a k-sparse unit vector x̂ achieving constant squared correlation5 with the signal x as soon as λ & √ k · log(np/k). That is, whenever λ & max‖x‖=1,‖x‖0=k〈W ,x ⊗p〉. Unfortunately, this approach runs in time exponential in k and takes super-polynomial time when p . k.6 As such, we assume p ≤ k from now on. Taking into account computational aspects, the picture changes. A good starting point to draw intuition for SSTM is the literature on sparse PCA and tensor PCA. We briefly outline some known results here. To simplify the discussion, we hide absolute constant multiplicative factors using O(·), Ω(·), ., and &, and hide multiplicative factors logarithmic in n using Õ(·).

heading: 1.1.1 Sparse PCA (Wigner noise)
text: Sparse PCA with Wigner noise exhibits a sharp phase transition in the top eigenvalue of Y for λ ≥ √ n [FP07]. In this strong signal regime, the top eigenvector7 v of Y correlates8 with x with high probability, thus the following spectral method achieves the same guarantees as the exhaustive search suggested above: compute a leading eigenvector of Y and restrict it to the top k largest entries in absolute value. Conversely, when λ < √ n, the top eigenvector of Y does not correlate with the signal x. In this weak signal regime, [JL09] proposed a simple algorithm known as diagonal thresholding: compute the top eigenvector of the principal submatrix defined by the k largest diagonal entries of Y . This algorithm recovers the sparse direction when λ & Õ(k), thus requiring almost an additional √ k factor when compared to inefficient algorithms. More refined polynomial-time algorithms (low-degree polynomials [dKNS20], covariance thresholding [DM16] and the basic SDP relaxation [dGJL07, dKNS20]) only improve over diagonal thresholding by a logarithmic factor in the regime n1−o(1) . k2 . n. Interestingly, multiple results suggest that this information-computation gap is inherent to the sparse PCA problem [BR13a, BR13b, DKWB19, dKNS20]. Subexponential time algorithms and lower bounds have also been shown. For instance, [DKWB19, HSV20] presented smooth trade-offs between signal strength and running time.9
5One could also aim to find a unit vector with correlation approaching one or, in the restricted setting of x ∈ {0,±1/ √ k}, aim to recover the support of x. At the coarseness of our discussion here, these goals could be considered mostly equivalent. 6Note that the problem input is of size np. So when p & k, exhaustive search takes nO(p) time which is polynomial in np. Thus, the interesting parameter regimes occur when p . k. 7By “top eigenvector” or “leading eigenvector”, we mean the eigenvector corresponding to “largest (in absolute value) eigenvalue”. 8More precisely, the vector consisting of the k largest (in absolute value) entries of v. 9Both works studied the single spike matrix setting. [HSV20] only considers the Wishart noise model and thus its guarantees cannot be compared to ours. [DKWB19] studied both the Wishart and Wigner noise models. In the Wishart noise model setting, both [HSV20] and [DKWB19] observe the same tradeoff between running time and signal-to-noise ratio. In the Wigner noise model setting, our algorithm and the algorithm of [DKWB19] offer the same smooth-trade off between running time and signal strength, up to universal constants.

heading: 1.1.2 Tensor PCA
text: In tensor settings, computing max‖x‖=1〈Y , x⊗p〉 is NP-hard already for p = 3 [HL13]. For even tensor powers p, one can unfold the tensor Y into a np/2-by-np/2 matrix and solve for the top eigenvector [MR14]. However, this approach is sub-optimal for odd tensor powers. For general tensor powers p, a successful strategy to tackle tensor PCA has been the use of semidefinite programming [HSS15, BGL16, HSS19]. Spectral algorithms inspired by the insight of these convex relaxations have also been successfully applied to the problem [SS17]. These methods succeed in recovering the single-spike x when λ & Õ ( np/4 ) , thus exhibiting a large gap when compared to exhaustive search algorithms. Matching lower bounds have been shown for constant degrees in the Sum-of-Squares hierarchy [BGL16, HKP+17] and through average case reductions [BB20].

heading: 1.1.3 Sparsity-exploiting algorithms and tensor algorithms
text: It is natural to ask how do the characteristics of sparse PCA and tensor PCA extend to the more general setting of SSTM. In particular, there are two main observations to be made. The first observation concerns the sharp computational transition that we see for k . √ n in sparse PCA. In these highly sparse settings, the top eigenvector of Y does not correlate with the signal x and so algorithms primarily based on spectral methods fail to recover it. Indeed, the best known guarantees are achieved through algorithms that crucially exploit the sparsity of the hidden signal. These algorithms require the signal strength to satisfy λ ≥ Õ( √ k), with only logarithmic dependency on the ambient dimension. To exemplify this to an extreme, notice how the following algorithm can recover the support of xxT with the same guarantees as diagonal thresholding, essentially disregarding the matrix structure of the data: zero all but the k2 largest (in absolute value) entries of Y . A natural question to ask is whether a similar phenomenon may happen for higher order tensors. In the highly sparse settings where k . √ n, can we obtain better algorithms exploiting the sparsity of the hidden vector? Recently, a partial answer appeared in [LZ20] with a polynomial time algorithm recovering the hidden signal for λ ≥ Õ(p · kp/2), albeit with suboptimal dependency on the tensor order p. The second observation concerns the computational-statistical gap in the spiked tensor model. As p grows, the gap between efficient algorithms and exhaustive search widens with the polynomial time algorithms requiring signal strength λ ≥ Õ ( np/4 ) while exhaustive search succeeds when λ ≥ Õ( √ n) [MR14]. The question here is: how strong is the dependency on p for efficient algorithms in sparse signal settings?
In this work, we investigate these questions in the high order tensors regime p ∈ ω(1). We present a family of algorithms with a smooth trade-off between running time and signal-to-noise ratio. Even restricting to polynomial-time settings, our algorithms improve over previous results. Furthermore, through the lens of low-degree polynomials, we provide rigorous evidence of an exponential gap in the tensor order p between algorithms and lower bounds. Remark. The planted sparse densest sub-hypergraph model [CPMB19, BCPS20, CPSB20] is closely related to SSTM. We discuss this model in Appendix E.

heading: 1.2 Results

heading: 1.2.1 Single spike setting
text: Consider first the restricted, but representative, case where the planted signal is a (k,A)-sparse unit vector with k non-zero entries having magnitudes in the range [ 1
A √ k , A√ k
] for some constant A ≥ 1.
We say that the signal is flat when A = 1 and approximately flat when A ≥ 1. Our first result is a limited brute force algorithm – informally, an algorithm that smoothly interpolates between some brute force approach and some “simple” polynomial time algorithm – that exactly recovers the signal support of the planted signal10. Theorem 1 (Algorithm for single spike sparse tensor PCA, Informal). Let A ≥ 1 be a constant. Consider the observation tensor
Y = W + λx⊗p
10A similar algorithm was analyzed by [DKWB19] for the special case of p = 2 and r = 1.
where the additive noise tensorW ∈ ⊗pRn contains i.i.d. N(0, 1) entries and the signal x ∈ Rn is a (k,A)-sparse unit vector with signal strength λ > 0. Let 1 ≤ t ≤ k be an integer. Suppose that
λ & √ t ( 2A2k
t
)p lnn .
Then, there exists an algorithm that runs in O(pnp+t) time and, with probability 0.99, outputs the support of x.
Let’s first consider Theorem 1 in its simplest setting where A = 1 and t is a fixed constant. For k . √ n, the theorem succeed when λ ≥ Õ(kp/2), thus improving over the guarantees of known
tensor PCA methods which require λ ≥ Õ(np/4). In addition, since support recovery is exact, one can obtain a good estimate11 of the planted signal by running any known tensor PCA algorithm on the subtensor corresponding to its support. Indeed, the resulting subtensor will be of significantly smaller dimension and the requirement needed on the signal strength by single-spike tensor PCA algorithms are weaker than the requirement we impose on λ (see Remark 9 for details). As a result, our algorithm recovers the guarantees of diagonal thresholding in the matrix (p = 2) setting. Our polynomial-time algorithm also improves over the result of [LZ20], which required λ & √ pkp log n, by removing the polynomial dependency of the tensor order p in the signal strength λ.12
Consider now the limited brute force parameter t. From the introductory exposition, we know that one can obtain a statistically optimal algorithm by performing a brute force search over the space of k-sparse flat vectors in Rn. The limited brute force algorithm is a natural extension that takes into account computational constraints by searching over the smaller set of t-sparse flat vectors, for 1 ≤ t ≤ k, to maximize 〈Y , u⊗p〉. The parametric nature of the algorithm captures both the brute force search algorithm (when t = k) and the idea of diagonal thresholding (when t = 1 and p = 2). As long as t ≤ k, using a larger t represents a direct trade-off between running time and the signal-to-noise ratio. Extending the result to approximately flat vectors, the dependency on A in the term ( 2A2 )p can be removed by increasing the computational budget to some value t′ ≥ 2A2t.

heading: 1.2.2 Multiple spikes
text: Theorem 2 (Algorithm for multi-spike sparse tensor PCA, Informal). Let A ≥ 1 be a constant. Consider the observation tensor
Y = W + r∑ q=1 λqx ⊗p (q)
where the additive noise tensor W ∈ ⊗pRn contains i.i.d. N(0, 1) entries and the signals x(1), . . . , x(r) ∈ Rn are (k,A)-sparse unit vectors with disjoint supports and corresponding signal strengths λ1 ≥ . . . ≥ λr > 0. Let 1 ≤ t ≤ k be an integer and 0 < ≤ 1/2. Suppose that
λr & 1 · √ t ( 2A2k
t
)p lnn and λr & A2p · (2 )p−1 · λ1 .
Then, there exists an algorithm that runs in O(rpnp+t) time and, with probability 0.99, outputs the individual signal supports of x(π(1)), . . . , x(π(r)) for some unknown bijection π : [r]→ [r].
Theorem 2 requires two assumptions on the signals: (1) signals have disjoint support; (2) there is a bounded signal strength gap of λr & A2p · (2 )p−1 · λ1. In the context of sparse PCA, algorithms that recover multiple spikes (e.g. [JL09, DM16]) only require the sparse vectors to be orthogonal. However, their guarantees are of the form λr ≥ Õ
(∣∣∣⋃q∈[r] supp (x(q))∣∣∣). That is, when the r signals have disjoint supports, they require the smallest signal to satisfy λr ≥ Õ (k · r). In comparison, already for constant t, Theorem 2 successfully recovers the supports when λr ≥ Õ(k),
11Recovery is up to a global sign flip since 〈u, v〉p = 〈u,−v〉p for even tensor powers p. 12The result of [LZ20] extends to the settings where Y = W + λX for an approximately flat tensor
X ∈ ⊗pRn. Both Theorem 1 and Theorem 2 can also be extended to these settings (see Appendix B.2).
thus removing the dependency on the number of signals and improving the bound by a 1/r factor13. Meanwhile, the bounded signal strength gap assumption is a common identifiability assumption (e.g. see [CMW13, DM16]). We remark that Theorem 2 provides a tradeoff between this signal strength gap assumption and the signal strengths: we can recover the supports with a smaller gap if the signal strengths are increased proportionally – increasing λr by a multiplicative factor α enables the algorithm to succeed with gap that is smaller by a multiplicative factor of 1/α. As an immediate consequence, we also obtain a tradeoff between gap assumption and running time: every time we double t (while ensuring 1 ≤ t ≤ k), λr increases by a factor of (1/ √ 2)p−1 and thus the algorithm can succeed with a smaller gap. Finally, as in the single spike case, the exact support recovery allow us to obtain good estimate of each signal by running known tensor PCA algorithms.
Remark We remark that these results can be extended to the general tensor settings
Y = W + r∑ q=1 λqX(q)
where for q ∈ [r], X(q) = x(q,1) ⊗ · · · ⊗ x(q,p) ∈ ⊗pRn in a natural way. See Appendix B.2.

heading: 1.2.3 An exponential gap between lower bounds and algorithms
text: SSTM generalizes both sparse PCA and tensor PCA. Hence, a tight hardness result for the model is interesting as it may combine and generalize the known bounds for these special cases. Here, we give a lower bound for the restricted computational model captured by low-degree polynomials. Originally developed in the context of the sum of squares hierarchy, this computational model appears to accurately predict the current best-known guarantees for problems such as sparse PCA, tensor PCA, community detection, and planted clique (e.g. see [HS17, HKP+17, Hop18, BHK+19, DKWB19, KWB19, dKNS20]).
Theorem 3 (Lower bound for low-degree polynomials, Informal). Let 1 ≤ D ≤ 2n/p and ν be the distribution of Z ∈ ⊗pRn with i.i.d. entries from N(0, 1). Then, there exists a distribution µ over tensors Y ∈ ⊗pRn of the form
Y = W + λx⊗p
where W ∈ ⊗pRn is a noise tensor with i.i.d. N(0, 1) entries, the marginal distribution of x is supported on vectors with entries { ±1/ √ k, 0 }n
, and x and W are distributionally independent, such that whenever
λ .
√ D
2p min
{( n
pD
)p/4 , ( k
pD
( 1 + ∣∣∣∣ln(npDek2 )∣∣∣∣))p/2 } ,
µ is indistinguishable14 from ν with respect to all polynomials of degree at most D.
Theorem 3 states that for certain values of λ, low-degree polynomials cannot be used to distinguish between the distribution of Y andW as typical values of low-degree polynomials are the same (up to a vanishing difference) under both distributions. The theorem captures known results in both sparse and tensor PCA settings. When p = 2, our bound reduces to λ . min {√ n, k√
D
( 1 + ∣∣ln ( 2nDek2 )∣∣)}, matching known low-degree bounds of [DKWB19] in the sparse PCA setting. Meanwhile, in the
tensor PCA settings (p ≥ 2, k = n), Theorem 3 implies a bound of the form λ . √ D ( n pD )p/4 , thus recovering the results of [KWB19].
13It is an intriguing question whether an improvement of 1/r can be achieved in the more general settings of orthogonal spikes. Our approach relies on the signals having disjoint support and we expect it to not be generalizable to orthogonal signals. This can be noticed in the simplest settings with brute-force parameter t = 1 and p = 2 where the criteria of Algorithm 3 for finding an entry of a signal vector is to look at the diagonal entries of the data matrix. In this case, the algorithm may be fooled since the largest diagonal entry can depend on more than one spike. Nevertheless, we are unaware of any fundamental barrier suggesting that such guarantees are computationally hard to achieve.
14In the sense that for any low-degree polynomial p(Y ) we have Eµp(Y )−Eνp(Y )√ Vp(Y ) ∈ o(1). See Appendix A.4.2.
For constant power p and k . √ n, our lower bound suggests that no estimator captured by polynomials of degree D . t log n can improve over our algorithmic guarantees by more than a logarithmic factor. However, for p ∈ ω(1), an exponential gap appears between the bounds of Theorem 3 and state-of-the-art algorithms (both in the sparse settings as well as in the dense settings).15 As a concrete example, let us consider the setting where p = n0.1 < k. The polynomial time algorithm of Theorem 1 requires λ ≥ Õ(kp/2) while according to Theorem 3 it may be enough to have λ ≥ Õ ( k/n0.1 )p/2 . Similarly, for k & √ np, known tensor algorithms recovers the signal for
λ ≥ Õ(np/4) while our lower bound only rules out algorithms for λ ≤ Õ ( n0.9·p/4 ) .
Surprisingly, for the distinguishing problem considered in Theorem 3, these bounds appear to be tight. For a wide range of parameters (in both the dense and sparse settings) there exists polynomial time algorithms that can distinguish the distributions ν and µ right at the threshold considered in Theorem 3 (see Appendix C). It remains a fascinating open question whether sharper recovering algorithms can be designed or stronger lower bounds are required.
Finally, we would like to highlight that this non-trivial dependency on p is a purely computational phenomenon as it does not appear in information-theoretic bounds (see Appendix D).
Remark Note that Theorem 3 is not in itself a lower bound for the recovery problem. However, any algorithm which obtains a good estimation of the signal vector x for signal strength λ ≥ √ k log n can be used to design a probabilistic algorithm which solve the distinguishing problem for signal strength Op(λ). Let us elaborate. Consider an algorithm that given Y = W + λx⊗p outputs a vector x̂ such that |〈x̂, x〉| ≥ 0.9. With high probability, max|z|2=1,|z|0=k |〈W , z ⊗p〉| ≤ Õ( √ k) and thus
|〈Y , x̂⊗p〉| ≥ λ · (0.9)p − Õ( √ k). Therefore, one can solve the distinguishing problem as follows: output “planted” if |〈Y , x̂⊗p〉| & √ k log n and “null” otherwise.

heading: 1.3 Notation and outline of paper
text: We write random variables in boldface and the set {1, . . . , n} as [n]. We hide absolute constant multiplicative factors and multiplicative factors logarithmic in n using standard notations: O(·), Ω(·), ., &, and Õ(·). We denote by e1, . . . , en ∈ Rn the standard basis vectors. For x ∈ Rn, we use supp(x) ⊆ [n] to denote the set of support coordinates. We say that x is a (k,A)-sparse vector if k ∈ [n], constant A ≥ 1, |supp(x)| = k, and 1
A √ k ≤ |x`| ≤ A√k for ` ∈ supp(x). When A = 1, we
say that x is a k-sparse flat vector and may omit the parameter A. For general A ≥ 1, we say that x is approximately flat. For an integer t ≥ 1, we define Ut = { u ∈ { − 1√
t , 0, 1√ t
}n : |supp(u)| = t } as the set of t-sparse flat vectors. For a tensor T ∈ ⊗pRn and a vector u ∈ Rn, their inner product is defined as 〈T, u⊗p〉 =
∑ i1,...,ip∈[n] Ti1,...,ipui1 . . . uip .
The rest of the paper is organized as follows: In Section 2, we introduce the main ideas behind Theorem 1 and Theorem 2. In Section 3, we flesh out some concrete unresolved research questions. Appendix A contains preliminary notions. We formally prove Theorem 1 and Theorem 2 in Appendix B. The lower bound Theorem 3 is given in Appendix C. We present an information theoretic bound in Appendix D. Appendix E discusses the planted sparse densest sub-hypergraph model. Finally, Appendix F contains any deferred technical proofs required throughout the paper.

heading: 2 Recovering signal supports via limited brute force searches
text: We describe here the main ideas behind our limited brute force algorithm. We consider the model
Model 4 (Sparse spiked tensor model). For A ≥ 1, r ≥ 1, k ≤ n we observe a tensor of the form
Y = W + r∑ q=1 λqx ⊗p (q) ∈ ⊗ pRn
15In particular, in the sparse settings k ≤ √np, the p−p/2 factor could not be seen in the restricted case of sparse PCA (as this factor is a constant when p = 2).
where W ∈ ⊗pRn is a noise tensor with i.i.d. N(0, 1) entries, λ1 ≥ . . . ≥ λr > 0 are the signal strengths, and x(1), . . . , x(r) are k-sparse flat unit length signal vectors.
We first look at the simplest setting of a single flat signal (i.e. A = 1 and r = 1). Then, we explain how to extend the analysis to multiple flat signals. For a cleaner discussion, we assume here that all the non-zero entries of the sparse vector x and vectors in the set Ut have positive sign. Our techniques also easily extend to approximately flat vectors (where A ≥ 1) and general signal tensors x(1) ⊗ · · · ⊗ x(p) ∈ ⊗pRn. We provide details for these extensions in Appendix B.

heading: 2.1 Single flat signal
text: Limited-brute force As already mentioned in the introduction, a brute force search over Uk for the vector maximizing 〈Y , u⊗p〉 returns the signal vector x (up to a global sign flip) with high probability whenever λ & √ k log n. This algorithm provides provably optimal guarantees but requires exponential time (see Appendix D for an information-theoretic lower bound). The idea of a limited brute force search is to search over a smaller set Ut (1 ≤ t ≤ k) instead, and use the maximizer v∗ to determine the signal support supp (x). The hope is that for a sufficiently large signal-to-noise ratio, this t-sparse vector v∗ will still be non-trivially correlated with the hidden vector x. Indeed as t grows, the requirement on λ decreases towards the information-theoretic bound, at the expense of increased running time.
As a concrete example, consider the matrix settings (p = 2). It is easy to generalize the classic diagonal thresholding algorithm ([JL09]) into a limited brute-force algorithm. Recall that diagonal thresholding identifies the support of x by picking the indices of the largest k diagonal entries of Y . In other words, the algorithm simply computes 〈Y , e⊗2i 〉 for all i ∈ [n] and returns the largest k indices. From this perspective, the algorithm can be naturally extended to t > 1 by computing the( k t ) vectors u ∈ Ut maximizing 〈Y , u⊗2〉 and reconstructing the signal from them. For t = k, the algorithm corresponds to exhaustive search.
With this intuition in mind, we now introduce our family of algorithms, heavily inspired16 by [DKWB19]. We first apply a preprocessing step to obtain two independent copies of the data.
Algorithm 1 Preprocessing Input: Y . Sample a Gaussian tensor Z ∈ ⊗pRn where each entry is an i.i.d. standard Gaussian N(0, 1). Return two independent copies Y (1) and Y (2) of Y as follows:
Y (1) = 1√ 2 (Y +Z) and Y (2) = 1√ 2 (Y −Z)
Algorithm 1 effectively creates two independent copies of the observation tensor Y . To handle the noise variance, the signal-to-noise ratio is only decreased by the constant factor 1/ √ 2. For simplicity, we will ignore this constant factor in the remainder of the section.
Algorithm 2 Single spike limited brute force Input: k, t and Y (1),Y (2) obtained from Algorithm 1. Compute v∗ := argmaxu∈Ut〈Y
(1), u⊗p〉. Compute the vector α ∈ Rn with entries α` := 〈Y (2), v⊗p−1∗ ⊗ e`〉 for every ` ∈ [n]. Return the indices of the largest k entries of α.
16The algorithm in [DKWB19] is a specialization of ours (with comparable guarantees) in the simplest setting of p = 2 and a single spike. However, looking at [DKWB19], it is a priori unclear how to generalize the result to the settings of our interest. This is especially true in the tensor settings (p ≥ 3) with multiple spikes, where signals may interfere with each other.
The signal support recovery process outlined in Algorithm 2 has two phases. In the first phase, we search over Ut to obtain a vector v∗ that is correlated with the signal x. In the second phase, we use v∗ to identify supp(x). The correctness of the algorithm follows from these two claims:
(i) The t-sparse maximizer v∗ shares a large fraction of its support coordinates with signal x.
(ii) The k largest entries of α belong to the support supp(x) of signal x.
Crucial to our analysis is the following standard concentration bound on Gaussian tensors. We directly use Lemma 5 in our exposition here, and formally prove a more general form in Appendix F.1.
Lemma 5. Let p ≤ n, t > 0 be an integer, andW ∈ ⊗pRn be a tensor with i.i.d. N(0, 1) entries. Then, with high probability, for any u ∈ Ut,
〈W ,u⊗p〉 . √ t log n .
For some constant 0 < ≤ 1/2, suppose that
λ & 1
· (1− )p−1 ·
√ t ( k
t
)p log n . (1)
For any u ∈ Ut with support supp(u) ⊆ supp(x), we have
〈Y (1), u⊗p〉 = λ〈x, u〉p + 〈W (1), u⊗p〉 ≥ λ · ( t
k
) p 2 −O (√ t log n ) .
On the other hand, any u ∈ Ut with support satisfying |supp(u) ∩ supp(x)| ≤ (1− ) · t has small correlation with Y (1) in the sense that
〈Y (1), u⊗p〉 = λ〈x, u〉p + 〈W (1), u⊗p〉 ≤ λ · (1− )p · ( t
k
) p 2 +O (√ t log n ) .
By Eq. (1), with high probability, v∗ will have at least a fraction (1− ) of its support contained in supp(x), yielding the first claim. Observe that v∗ does not completely overlap with x. A priori, this might seem to be an issue. However, it turns out that we can still use v∗ to exactly reconstruct the support of x. Indeed, for all ` ∈ supp(x),
α` = λ · x` · 〈x,v∗〉p−1 + 〈W (2), v⊗p−1∗ ⊗ e`〉
≥ λ · (1− ) p−1
√ k
· ( t
k
) p−1 2
+ 〈W (2), v⊗p−1∗ ⊗ e`〉
& 1 · √
log n+ 〈W (2), v⊗p−1∗ ⊗ e`〉 .
Now, by independence ofW (2) and v∗, 〈W (2), v⊗p−1∗ ⊗ e`〉 behaves like a standard Gaussian. Thus, with high probability, ∣∣∣〈W (2), v⊗p−1∗ ⊗ e`〉∣∣∣ . √log n and α` & √log n. Conversely, if ` is not in the support of the signal, then α` . √ log n. So, the vector α acts as indicator of the support of x!
Remark 6. In its simplest form of t = 1, Algorithm 2 does not exploit the tensor structure of the data: it performs entry-wise search for the largest (in magnitude) over a subset of Y . However, this is no longer true as t grows. For t = k, the algorithm computes the k-sparse flat unit vector u maximizing 〈Y (1), u⊗p〉.

heading: 2.2 Multiple flat signals with disjoint signal supports
text: Consider now the setting with r > 1 spikes. Recall that we assumed the vectors x(1), . . . , x(r) to have non-intersecting supports. We also assumed that for any q, q′ ∈ [r] and some fixed scalar 0 ≤ κ ≤ 1, if λq ≥ λq′ , then λq′ ≥ κ · λq. We remark that we may not recover the signal supports in a known
order, but we are guaranteed to recover all of them exactly. For simplicity of discussion, let us assume here that we recover the vector x(i) at iteration i.
The idea to recover the r spikes is essentially to run Algorithm 2 r times. At first, we compute the t-sparse vector v∗ by maximizing the product 〈Y (1), v⊗p∗ 〉. Then, using v∗, we compute the vector α to obtain a set I1 ⊆ [n]. With high probability, we will have I1 = supp ( x(1) ) and so we will exactly recover the support of x(1). In the second iteration of the loop, we repeat the same procedure with the additional constraint of searching only over the n− k dimensional subset of Ut containing vectors with disjoint support from I1. Similarly, at iteration i, we search over the subset of Ut containing vectors with disjoint support from
⋃ 1≤j<i Ij . As before, we first preprocess the data to create two
independent copies Y (1) and Y (2). Concretely:
Algorithm 3 Multi-spike limited brute force Input: k, t, r and Y (1),Y (2) obtained from Algorithm 1. Repeat for i = 1 to r:
Compute v∗ := argmaxu∈Ut〈Y (1), u⊗p〉 subject to supp (v∗) ∩ ( ⋃ 1≤j<i Ij ) = ∅. Compute the vector α ∈ Rn with entries α` := 〈Y (2), v⊗p−1∗ ⊗ e`〉 for every ` ∈ [n]. Let Ii be the set of indices of the largest k entries of α.
Return I1, . . . , Ir.
The proof structure is similar to that of Algorithm 2 and essentially amounts to showing that the claims (i) and (ii) described in Section 2.1 hold in each iteration.
Let λmin = minq∈[r] λq and λmax = maxq∈[r] λq. For some 0 < ≤ 1/2, let κ & (
1−
)p−1 such
that λmin ≥ κ · λmax. Suppose that
λmin & 1
· (1− )p ·
√ t ( k
t
)p log n and λmin & (
1−
)p−1 · λmax . (2)
Consider an arbitrary iteration i and suppose that we exactly recovered the support of one signal in each of the previous iterations. Without loss of generality, assume that λmax is the largest signal strength among the yet to be recovered signals, and let x(max) be one such corresponding signal.
For u ∈ Ut satisfying supp(u) ⊆ supp(x(max)), we have
〈Y (1), u⊗p〉 = λmax〈x(max), u〉p + 〈W (1), u⊗p〉 ≥ λmax · ( t
k
) p 2 −O (√ t log n ) .
On the other hand, for any u ∈ Ut such that ∣∣supp (u) ∩ supp (x(q))∣∣ ≤ (1− ) · t for all q ∈ [r],
〈Y (1), u⊗p〉 = ∑ q∈[r] λq〈x(q), u〉p + 〈W (1), u⊗p〉
≤ λmax · ( t
k
) p 2 · ((1− )p + p) +O (√ t log n )
≤ λmax · ( t
k
) p 2 · (1− )p−1 +O (√ t log n ) .
Thus, as in Section 2.1, it follows that v∗ satisfies ∣∣supp (v∗) ∩ supp (x(i))∣∣ ≥ (1− ) · t for some signal x(i). Note that x(i) may not be x(max). Even though v∗ does not exactly overlap with any of the signal vectors, we will not accumulate an error at each iteration. This is because, analogous to the single spike setting, we can exactly identify the support of a signal through α. For any
` ∈ supp ( x(i) ) , it holds that α` & √ log n as before because ∣∣∣〈W (2), v⊗p−1∗ ⊗ e`〉∣∣∣ . √log n. Conversely, since signal supports are disjoint, we see that for ` /∈ supp ( x(i) ) ,
α` = ∑ q∈[r] λq · x(q),` · 〈x(q), v∗〉p−1 + 〈W (2), v⊗p−1∗ ⊗ e`〉
≤ λmax · p−1√ k · ( t k ) p−1 2 +O (√ log n )
≤ λmin κ · p−1 √ k · ( t k ) p−1 2 +O (√ log n )
. √ log n .
So, once again, α exactly identifies the support of x(i) with high probability. Remark 7 (On the strength of the assumption on κ). As already briefly discussed in Section 1.2, the algorithm provides a three-way trade-off between signal gap κ, signal-to-noise ratio λ and running time. By appropriately choosing the constant > 0, the algorithm can tolerate different values of
κ. Indeed, the above analysis holds as long as κ & (
1−
)p−1 . This suggests two ways in which
we can loosen the requirement λmin ≥ κ · λmax and still successfully recover the signals through Algorithm 3. One is increase the running time, so that we can decrease without increasing the signal-to-noise ratio λmin. The other is to decrease and increase the value of λmin accordingly. Remark 8 (On independent copies of Y ). To clarify why it suffices to have 2 independent copies of Y even for multiple iterations, observe that at each iteration i, the choice of the set Ii depends only on the vector v∗ with high probability. Consider the following thought experiment where we are given a fresh copy Y (i) of Y in the second phase of each iteration i of the algorithm (while still using only a single copy Y (1)for all the first phases). Even with fresh randomness, the result is the same as Algorithm 3 with high probability because at each iteration the choice of maximizer v∗ causes the same output. Remark 9 (Reconstructing the signals from their supports). After recovering individual signal supports, one can reconstruct signals using known tensor PCA algorithms (e.g. [MR14, HSS15]) on the subtensor defined by each recovered support. The signal strength required for this new subproblem is weaker and is satisfied by our recovery assumptions. For instance, by concatenating our algorithm with [HSS15, Theorem 7.1], one obtains vectors x̂(1), . . . , x̂(r) such that
∣∣〈x̂(i), x(i)〉∣∣ ≥ 0.99, for any i ∈ [r], with probability 0.99.

heading: 3 Open questions
text: Open question 1. Theorem 2 improves over existing sparse PCA multi-spike recovery algorithms (which only assumed orthogonal spikes) by a factor of 1/r in the case where these planted signals have disjoint support. Can one still obtain an improvement of 1/r if we only assume orthogonality? Open question 2. For p ∈ ω(1), there is an exponential gap between the bounds of Theorem 3 and state-of-the-art algorithms. A natural question is whether one can design better recover algorithms or prove stronger lower bounds for this range of tensor power p ∈ ω(1)?
Remark (Societal impact). This work does not present any foreseeable negative societal consequence.
Acknowledgments and Disclosure of Funding
This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 815464). This research/project is supported by the National Research Foundation, Singapore under its AI Singapore Programme (AISG Award No: AISG-PhD/2021-08-013). We thank David Steurer for several helpful conversations. We thank Luca Corinzia and Paolo Penna for useful discussions about the planted sparse densest sub-hypergraph model.

Review summaries:
This paper studies the problem of sparse tensor PCA (a generalization of sparse PCA and tensor PCA). The main result is an algorithm that -- for the highly sparse regime in which the sparsity is at most the square root of the signal dimension -- naturally interpolates between a poly-time efficient algorithm and exponential exhaustive search. This result is complemented by lower bounds for the rather popular computational model captured by low-degree polynomials. The lower bounds match the algorithm guarantees (although for the different problem of distinguishability) in the regime of constant tensor order $p$.   The results of this paper recover several results of the literature, and provide a non-trivial extension. 

# Review for "The Complexity of Sparse Tensor PCA"  This paper studies the "Sparse Tensor PCA" problem, which is a common generalization of two well-studied problems: sparse principal component analysis and tensor principal component analysis. The problem is as follows: given a "single-spike" tensor of the form  $T = \lambda \cdot x^{\otimes p} + W$  where $W$ has i.i.d. Gaussian entries, $\lambda > 0$ is the "signal strength", and $x \in R^n$ is $k$-sparse for some $k \ll n$, find $x$. (The problem can be generalized to more than one spike $x$; this is treated carefully in the paper.)  The paper contains interesting new algorithms and computational lower bounds for this problem. 

The authors propose a family of algorithms for solving the sparse tensor PCA problems, which is based on smooth interpolation between a polynomial-time algorithm and exponential-time search algorithm. Improved results are achieve for sparce PCA in the lower signal/noise ratio regime in distinct sparse signal cases, compared with existing guarantees. Lower bound analysis is also provided for the proposed algorithm.

This paper focuses on the tradeoff between computational complexity and conditions for recovery in a symmetric sparse tensor estimation problem. There has been a great deal of work on related problems recently in the case of matrices (order 2 tensors) and the case of tensors (without necessarily imposing sparsity constraints). The main achievability result (Theorem 1) provides sufficient condition for single spike model in terms of the ambient dimension, the order the tensors, the degree of sparsity, the signal strength, and an integer parameter the interpolates between a ````"low complexity" algorithm  "high complexity" algorithm (essentially brute force search). Similar results of for a multiple spike model are given in Theorem 2.   Necessary conditions for a certain testing problem are considered in Theorem 3 which provides necessary conditions such that a polynomial of bounded degree can "distinguish" between a planted model (with a signal of a given strength) and a null mull of only Gaussian noise. Here the criterion for distinguishability is a certain L2 test.  While similar results have been obtained in the literature, the authors argue that their results improve upon the existing ones in a number of ways. 

This paper studies the estimation of a rank-one (or more generally a low-rank) tensor $x^{\otimes p}$ ($x \in \mathbb{R}^n$ is assumed to have unit norm) given $Y = \lambda x^{\otimes p} + W$, where $\lambda > 0$ and $W$ is a Gaussian noise tensor. Here the signal vector $x$ is assumed to have only $k$ non-zero entries.  The main contribution of the paper is an algorithm that recovers with high probability the support of $x$ in time $O(n^{p+t})$ provided that $\lambda \geq O(\sqrt{t} (k/t)^{p/2})$. The number $t$ can be arbitrarily chosen between $1$ and $k$, illustrating a tradeoff between statistical constraints on the signal-to-noise ratio $\lambda$ and runtime constraints on the estimation algorithm. This generalizes to the 'tensor settings' previous results obtained about sparse PCA.  The authors also prove an information-theoretic lower bound for low-degree polynomials. More precisely, they prove that no test based on polynomial of degree at most $D$ can distinguish between the planted model $Y = \lambda x^{\otimes p} + W$ and the pure noise model $Y=W$, when $\lambda$ is smaller than some quantity depending on $D, n, k, p$. This bound generalizes previously obtained results for PCA and tensor PCA. This bound is tight for low-degree polynomials. However there is a gap with the  bound $\lambda \geq O(\sqrt{t} (k/t)^{p/2})$ needed for the algorithm that recovers the support.